pyGIMLi: An open-source library for modelling and inversion in geophysics

Many tasks in applied geosciences cannot be solved by single measurements, but require the integration of geophysical, geotechnical and hydrological methods. Numerical simulation techniques are essential both for planning and interpretation, as well as for the process understanding of modern geophysical methods. These trends encourage open, simple, and modern software architectures aiming at a uniform interface for interdisciplinary and ﬂ exible modelling and inversion approaches. We present pyGIMLi (Python Library for Inversion and Modelling in Geophysics), an open-source framework that provides tools for modelling and inversion of various geophysical but also hydrological methods. The modelling component supplies discretization management and the numerical basis for ﬁ nite-element and ﬁ nite-volume solvers in 1D, 2D and 3D on arbitrarily structured meshes. The generalized inversion framework solves the minimization problem with a Gauss-Newton algorithm for any physical forward operator and provides opportunities for uncertainty and resolution analyses. More general requirements, such as ﬂ exible regularization strategies, time-lapse processing and different sorts of coupling individual methods are provided independently of the actual methods used. The usage of pyGIMLi is ﬁ rst demonstrated by solving the steady-state heat equation, followed by a demonstration of more complex capabilities for the combination of different geophysical data sets. A fully coupled hydrogeophysical inversion of electrical resistivity tomography (ERT) data of a simulated tracer experiment is presented that allows to directly reconstruct the underlying hydraulic conductivity distribution of the aquifer. Another example demonstrates the improvement of jointly inverting ERT and ultrasonic data with respect to saturation by a new approach that incorporates petrophysical relations in the inversion. Potential applications of the presented framework are manifold and include time-lapse, constrained, joint, and coupled inversions of various geophysical and hydrological data sets.


Introduction
In modern geophysical applications, it is often desired to maximize information on the subsurface by a combination of different measurement methods. When dynamic changes are monitored, an additional link to corresponding process models (e.g., hydrogeological or geomechanical models) can lead to an improved process understanding and offers opportunities to estimate multi-physical parameters of the subsurface. Joint and process-based inversions are therefore thriving research topics in the emerging field of hydrogeophysics (e.g., Binley et al., 2015;Linde and Doetsch, 2016).
However, such endeavors are associated with considerable technical challenges. The required coupling of different numerical models represents a potential impediment for many practitioners and students. Even technically versatile users commonly end up building individually tailored solutions by linking various existing (and potentially commercial) software packages through scripts, which hinders the reproducibility of scientific findings (Peng, 2011). This motivates and supports the need for open, simple, and modern software architectures for the main numerical tasks in geophysics. Uieda et al. (2013) present a software for geophysical data analysis with a focus on gravity and magnetic methods. Forward modelling routines are available to simulate gravitational and magnetic fields on various 2D and 3D meshes including tesseroids (spherical prisms). The inversion package enables non-linear parameter estimation with Levenberg-Marquardt, steepest-decent, as well as Gauss-Newton approaches. Hansen et al. (2013) provide a general inversion software for geophysical problems. Linear inverse Gaussian problems are solved using a least-squares solver, whereas general non-linear (i.e. non-Gaussian) inverse problems are solved with an extended Metropolis algorithm. While the software provides flexible inversion approaches for geophysical problems, its forward modelling functionality focuses on linear and non-linear travel time computations. Schaa et al. (2016) present a finite-element library for the solution of linear and non-linear, coupled, and time-dependent partial differential equations. The authors demonstrate the modelling capabilities based on 3D electrical resistivity and 2D magnetotelluric simulations. Cockett et al. (2015) present a software for simulation and gradient based parameter estimation in geophysics. Their approach is based on finite volume discretizations on structured and semi-structured meshes and includes convex optimization algorithms.
The authors state that joint and integrated inversions are generally possible by means of multiple misfit functions, physics, and regularization functionals.
The above-mentioned software packages are open-source, welldocumented, and aim at extensibility and reproducibility of geophysical simulations. With a comparable motivation in mind, we present pyGIMLi (Geophysical Inversion and Modelling Library in Python), a versatile and computationally efficient framework for modelling and inversion in geophysical applications. In contrast to existing approaches, pyGIMLi explicitly targets the modern aspirations in geophysics including constrained, joint, and process-based inversions together with the required forward modelling necessities.
pyGIMLi has been in active development since 2009 and offers modular functionality accessible from different levels of abstraction aiming at satisfying the diverse needs in research and education. The Python programming language was chosen as the basis for pyGIMLi for its free, flexible, and cross-platform-compatible nature, which therefore makes it widely used in the (geo)scientific community (e.g., Guyer et al., 2009;Logg and Wells, 2010;P erez et al., 2011;Wellmann et al., 2012;Uieda et al., 2013;Cockett et al., 2015;Weigand and Kemna, 2016;Hector and Hinderer, 2016;Schaa et al., 2016). One distinct advantage is that it can be easily extended by compiled modules from C or Fortran for example, allowing users to extend legacy code or outsource timeconsuming parts in computationally efficient extensions. We make use of this flexibility and have implemented all runtime sensitive parts in a Cþþ core library. Complete Python bindings to this core library are complemented by functionality written in pure Python, thus offering both efficiency and flexibility for the rapid development of robust modelling and inversion applications. The modelling component offers mesh management as well as finite-element and finite-volume solvers in 1D, 2D and 3D. The inversion component is based on a deterministic Gauss-Newton algorithm and works with any physical forward operator provided. Several post-processing routines are provided to visualize results in 2D using Matplotlib (Hunter, 2007) and in 3D using the software ParaView (Ayachit, 2015) or Mayavi (Ramachandran and Varoquaux, 2011).
After an introduction of the software architecture including detailed descriptions of the different abstraction levels, the generalized inversion framework is presented. This is followed by a demonstration on how to perform basic modelling. We then emphasize the main purpose of pyGIMLi, i.e. readily integrating interdisciplinary data sets, by two applications: i) a fully coupled hydrogeophysical inversion that directly inverts for the hydraulic conductivity distribution of the aquifer and ii) a new petrophysical joint inversion based on electrical resistivity and travel time tomography to directly estimate water saturation.

Software design
The main tasks to solve with pyGIMLi are modelling and inversion. The modelling component is realized with a finite element=volume toolbox with all common cell shapes for linear and quadratic base functions in 1D, 2D, and 3D domains. The default inversion is implemented with a generalized Gauss-Newton approach with flexible regularization.
In addition to the core library, the extension part of pyGIMLi also consists of third-party dependencies that provide advanced functionality. We apply mesh generators like Triangle (Shewchuk, 1996), Tetgen (Si, 2015) and Gmsh (Geuzaine and Remacle, 2009), as well as high-level numerical solvers like SuiteSparse (Davis, 2006). Fig. 1 shows the basic architecture of pyGIMLi with different abstraction levels written in Python. Based on the extension part, there are Python modules to exploit the easy-to-use scripting ability and to provide classes that are easier to maintain and to test. The equation level provides a general interface for the solution of common partial differential equations (PDE). The modelling level provides functionality of forward operators that aim for supporting specific geophysical methods.
In the application level, we define general frameworks to solve basic and advanced inversion tasks like time-lapse or joint inversion. The method managers combine all levels to solve the complete task for geophysical problems with specific data. All these abstraction levels communicate through unified interfaces and can thus be combined independently from the underlying geophysical method.

Cþþ core library
Python is a powerful scripting language for rapid development progress, but lacks runtime performance for pure Python code. It supports multiple paradigms that allow functional and object oriented programming and can be easily extended by pre-compiled high performance Cþþ extensions. We developed a Cþþ core library with a strong object oriented design for all runtime sensitive needs and provide full access to all parts of this core library by automatically generated Python bindings.
One main benefit of using Python is abstract prototyping, e.g., we implemented the main inverse solver in Cþþ with the use of an abstract base matrix and can use it with any custom advanced matrix directly from Python. This results in flexibility of the used matrix types with minimal coding effort and minimal runtime impairments. In addition to dense and different types of sparse matrices, there are specialized matrices like row or column-scaled matrices or Kronecker matrices. Furthermore, we implemented a block matrix containing references to a number of matrices of arbitrary type. It can be used for efficient constraint matrices or joint Jacobian matrices without performance loss.
For the high flexibility of regularization, we implemented a so called "region concept", where one can define parts of the inversion domain and treat, configure, or couple these regions individually. An overview of the background and different possibilities of incorporating information is given by Rücker (2011). Coscia et al. (2011) present an ERT inversion example, where different subsurface parts (geological layers, boreholes) have been decoupled and treated by different regularization operators.
The main performance drawback of the Python interpreted language is repeating code through nested loops due to the dynamic type conversion of the Python interpreter, which needs considerable time between two script instructions. A numerical package applying advanced geometry and mesh features spends significant runtime in mesh management, numerical integration, interpolation, or common tasks like nearest-neighbor search. Therefore, we implemented mesh management and related functionality as part of the Cþþ core library.
The majority of numerical approaches to the solution of partial differential equations (PDE) are based on a spatial discretization, called mesh or grid. A mesh combines N nodes, C cells, and B boundaries and represents the modelling domain Ω ¼ ∪ C i¼1 C i and its outer boundary A mesh can be imported from external mesh generators or generated by provided utility functions.
We consider a grid to be a structured form of a mesh (e.g., a regular discretization into quadrangles or hexahedrons) and treat them equally. Nodes N ¼ fN i ðrÞg with i ¼ 1…N represent N discrete position vectors r 2 ℝ 1 , ℝ 2 , or ℝ 3 . A cell C is a collection of nodes that spans the subdomain Ω C ¼ ∪N j in the same dimensional space of the nodes. As cell shapes we implemented simplexes such as triangles, quadrangles in ℝ 2 and tetrahedrons, hexahedrons, or prisms in ℝ 3 . Different cell types can be combined in one mesh, to combine structured discretization in the region of interest with progressively unstructured coarsening towards the boundaries for example. A boundary B is a collection of nodes spanning a subdomain Γ C ¼ ∂Ω C ¼ ∪N j representing the outer boundary of a cell. The boundary shapes arise as a result of the associated cell shape and are edges in ℝ 2 and triangles or quadrangles in ℝ 3 . For all shapes we provide numerical integration rules for several base functions to assemble the numerical prerequisites for finite element and finite volume analyses.

Equation level
The equation level provides an interface to solve common PDEs on a given mesh, which comprises all geometric specifications, e.g., topography or known subsurface structures. Currently the equation level provides entry points for solving the following two main PDE types, which cover a wide range of methods in applied geophysics from potential fields to wave propagation.
The finite element method (FEM) with linear or quadratic basis functions solves for u ¼ uðr; tÞ with r be the node positions by calling: u ¼ pygimli.solver.solveFiniteElement (mesh, a, b, f, uB, duB, t, u0) The finite volume method (FVM) with linear basis functions solves for u ¼ uðr; tÞ with r at the cell centers by calling: u ¼ pygimli.solver.solveFiniteVolume (mesh, a, b, v, f, uB, duB, t, u0) The solution space u ¼ uðr; tÞ can be 2 ℝjℂ for r 2 ℝ 1 ; ℝ 2 ; ℝ 3 . The coefficients a; b; c; v; u B ; u ∂B might be zero, scalars, arrays or functions(NjCjB, t) that are evaluated at runtime, e.g., they might be the solution of another solver run.
Note, that both solver functions are designed to give the user an easy as possible access to the modelling capabilities of the core library. However, most common elliptical, parabolic, and advection-type problems can be solved with these two functions, hyperbolic or curl-operator based systems can only be approached by nesting two parabolic equations. The accuracy of parabolic FVM is limited due to the linear base functions. Fortunately, highly specialized state-of-the art Python libraries exist (Guyer et al., 2009;Logg and Wells, 2010), which can be easily integrated if higher accuracy is needed.

Modelling level
The physics level represents a collection of classes to solve a forward or simulation task for a specific geophysical discipline by utilizing the equation level or applying suitable calculations. A forward operator (FOP) F ðmðr; tÞÞ maps a discrete parameter distribution m ¼ fm j g with j ¼ 1…M model parameters to the data vector d ¼ fd i g with i ¼ 1…N data. Each FOP is inherited from the base class ModellingBase implemented in the Cþþ core library and is accessed by the generalized interface: The method response needs to be filled with the appropriate numerical calculation. Note that the spatial or temporal discretization required in inversion is usually different to the numerical needs for accurate forward modelling. Therefore, pyGIMLi provides inter-and extrapolation tools to map parameters from one mesh to another or assigning them to specific regions.
The method createJacobian is the interface for calculating the entries of the Jacobian matrix: that might be needed for an inversion, but as it depends on the modelling problem, its part of the forward operator. It is beneficial to implement createJacobian with a method-specific Jacobian generation approach. However, if this function is not implemented, the base class ModellingBase provides a parallelized default mechanism to fill the entries for J using a finite-difference (brute-force) approach, i.e., repeated forward calculations with perturbed model parameters.
The implemented forward operators are summarized in Table 1:

Application level
The highest abstraction layer is the application level, and therefore it represents the first entrance point for scientific development and enduser interaction. It contains method-independent frameworks and method-specific managers with embedded apps.

Method managers and applications
The method managers represent the state of application that provides a full set of actions to proceed all tasks for a single discipline in applied geophysics and are considered entry points for end user interaction. These classes use a preconfigured instance of the generalized inversion framework and apply an appropriate forward operator with unified interfaces so that the methods simulate and invert are automatically available and controlled by keyword lists. Additionally, each method should provide the method-specific functions loadData, showData and showResults.
A typical use of a method manager to visualize and invert data for the example of ERT is: There are managers for practically all methods (Table 1) and combinations of them. Some are already including frameworks like LCI-type inversion (e.g., Costabel et al., 2016) or block-joint inversion (e.g., Günther and Müller-Petke, 2012). Moreover, there are managers involving simpler forward algorithms like fitting of complex conductivity/permittivity spectra (Loewer et al., 2016;Hupfer et al., 2016).

Inversion frameworks
Inversion frameworks are generalized, abstract approaches to solve a specific inversion problem without specifying the appropriate geophysical methods. This can be a specific regularization strategy, an alternative formulation of the inverse problem or algorithms of routine inversion. It is initialized by specific forward operators or managers that provide them. An example of how such an inversion framework is constructed for a petrophysical joint inversion is given in Appendix C.
2.4.2.1. Gauss-Newton inversion. The default inversion framework is based on the generalized Gauss-Newton method and is compatible with any given forward operator and thus applicable to various physical problems. We state the inversion problem as minimization of an objective function consisting of data misfit and model constraints: Note that we do not include inequality constraints in the minimization but use transformations to restrict parameters to reasonable ranges (e.g.,  Kim and Kim, 2011). W d is the data weighting matrix containing the inverse data errors, W m is the model constraint matrix (e.g., a first-order roughness operator), and m 0 is a reference model. The dimensionless factor λ scales the influence of the regularization term. There is a wide range of different regularization methods (different kinds of smoothness and damping, mixed operators, anisotropic smoothing) The existing ones can be used flexibly to constrain different model parameters or subsurface parts (regions), but also be extended by own functions. The application of the Gauss-Newton scheme on minimizing (4) yields the model update Δm k in the k th iteration (Park and Van, 1991): which is solved using a conjugate-gradient least-squares solver . The inversion process including the region-specific regularization is sketched in Fig. 2. All matrices of the inversion formulation can be directly accessed from Python and thereby offer opportunities for uncertainty and resolution analysis as well as experimental design (e.g., Wagner et al., 2015). Beyond different inversion approaches there are so-called frameworks for typical inversion (mostly regularization) tasks. Examples that are already implemented in pyGIMLi are for instance: Marquardt scheme inversion of few independent parameters, e.g., fitting of spectra (Loewer et al., 2016) Soil-physical model reduction incorporating soil-physical functions Costabel and Günther, 2014) Classical joint inversion of two data sets for the same parameter like DC and EM (Günther, 2013) Block joint inversion of several 1D data using common layers, e.g., MRSþVES (Günther and Müller-Petke, 2012) Sequential (constrained) inversion successive independent inversion of data sets, e.g., classic time-lapse inversion (e.g., Hübner et al., 2015) Simultaneous constrained inversion of data sets of data neighbored in space (LCI, e.g., Costabel et al., 2016), time (full time-lapse) or frequency  Structurally coupled cooperative inversion of disparate data based on structural similarity (e.g., Ronczka et al., 2017) Structure-based inversion using layered 2D models (Attwa et al., 2014)

Exemplary applications
To demonstrate the usage of pyGIMLi, several examples are given along with the code. For the sake of brevity, all examples use minimalist 2D geometries, but are directly transferable to 3D and more complex geometries, if a corresponding mesh is provided.

Simulating heat transfer based on a simple geometry
A simplistic modelling example shows the basic steps for solving a steady-state heat equation on the equation level. Fig. 3 shows the complete Python source code and the resulting images for model creation and finite element calculation. In the preamble (Fig. 3, lines 1,2) the necessary pyGIMLi namespace and the pyGIMLi mesh generation package are imported and abbreviated with the alias names pg and mt, respectively.
We assume three layers with a block inside the second layer. pyGIMLi provides basic geometry building utilities accessible through the package alias mt. The commands mt.createWorld and mt.createBlock create the desired geometric entities that are combined by the command mt.mergePLC (Fig. 3, line 12). The resulting geometry definition, so called piecewise linear complex (PLC), contains nodes, boundary elements, and region descriptions to represent the entire model geometry. Fig. 3 line: 13 creates an image of the model geometry using the pg.show command, which is the most straightforward way to view meshes, geometries, and data.
As we need a mesh, we forward the given PLC to the external mesh generator called Triangle (Shewchuk, 1996). The pg.show command is used again to view the resulting mesh ( Fig. 3 line: 15).
In the next step, we use solveFiniteElement from the equation level directly with the generated mesh to perform the FEM calculation. The arguments configure the requested PDE and control the underlying material parameters and boundary conditions. Most arguments are treated in a flexible manner, e.g., in this case a is a map that translates the four markers of the four geometry regions into a distribution of the thermal diffusivity a to create a layered background of a ¼ ½1; 2; 3 m 2 /s with the heterogeneous block (region 4) to be a ¼ 0:1 m 2 /s. The boundary conditions are controlled with the uB argument, e.g., as of Dirichlet type with a fixed temperature T ¼ 1 K for the boundary with marker 4 (bottom) and a fixed temperature T ¼ 0 K at the boundary with marker 3 (surface). The other boundaries obtain natural Neumann (noflow) boundary conditions by default.
The array for the resulting temperature distribution T can then be viewed by the calling pg.show again. The pg.show command can be also used to combine several plotting approaches to achieve the final result image of the given geometry with the estimated temperature distribution.

Fully coupled hydrogeophysical inversion
Monitoring of hydraulic processes with geophysical methods has become very popular. The classic way is to solve the geophysical problem, followed by later analysis to infer hydraulic properties. In recent years, a small number of approaches were presented that estimate hydraulic conductivity directly from geophysical observations (e.g., Pollock and Cirpka, 2010;Mboh et al., 2012;Camporese et al., 2015;Wagner, 2016). The small number of available approaches, however, mostly relies on customized and often proprietary software prohibiting reuse and advancement by other researchers. In the following, we show that coupled problems can be readily solved using pyGIMLi in a consistent and fully reproducible manner by a concatenation of hydraulic and geophysical simulations linked through petrophysical transformations.
We apply the geometry of the example from Sec. 3.1 and map hydraulic conductivities to the different regions of the model as shown in Fig. 4. The first layer of the model is considered a non-conducting topsoil followed by a conductive aquifer and a less conductive basement. The heterogeneous block inside the aquifer is assumed to represent a low conductive anomaly.
The fluid flow velocity in a porous medium of slow non-viscous and non-frictional hydraulic movement are governed by the Darcy equation (Whitaker, 1986): leading to with the hydraulic conductivity tensor K. We can solve equation (8) on the equation level for a given hydraulic potential p ¼ p 0 ¼ 0:75 m on the left and p ¼ 0 on the right boundary of the modelling domain, equaling a hydraulic gradient of 1.75%. The sought hydraulic velocity distribution can then be calculated as the gradient field of v ¼ À∇p. Extending the example with the following code fragment: leads to the velocity distribution displayed in Fig. 5.
The flow direction is from left to right and exhibits an increased velocity in the aquifer due to its larger hydraulic conductivity. The anomaly in the aquifer considerably interferes the velocity field and causes the field to circumvent this rectangular body.
In the next step we use this velocity field to simulate the dynamic movement of a particle (e.g., salt) concentration cðr; tÞ in the aquifer by using the advection-diffusion equation (e.g., Bechtold et al., 2012) as a result of a source S. The molecular diffusion coefficient D (in water ≈1⋅10 À9 m 2 /s) is negligible. However, in porous media a diffusion-like spreading characteristics called dispersion takes place that is governed by the same term. We choose a common velocity-depending dispersion coefficient D ¼ αjvj with a dispersivity α ¼ 1⋅10 À2 m (Bechtold et al., 2012). The particles are injected for six days at the position r s ðx ¼ À19:1; y ¼ À5:0Þ in the aquifer with an amount of S ¼ 0:001 g/ls. The whole simulation time is 12 days in summary and we used 1600 time steps to fulfill the Courant-Friedrichs-Lewy condition (Courant et al., 1967) ensuring the particle movement does not exceed cell dimensions within one time step. Solving equation (9) on the equation level with the finite volume solver results in a particle concentration cðr; tÞ (in g=l) for each cell center and time step. The following extension to the minimalist script: C. Rücker et al. Computers and Geosciences 109 (2017) 106-123 leads to the resulting temporal and spatial distribution of the concentration c and is shown in Fig. 6 after 0, 2, 4, 6, 8, and 10 days. We clearly see the particles spreading from the injection point and moving within the aquifer. The low conductive heterogeneous block disturbs the flow field and forces a large portion of the salt tracer to circumvent the obstacle. After 6 days the injection stops and the particle are washed out by the groundwater flow. The particles penetrate the heterogeneous block more slowly and remain longer due to the lower flow velocity.
If we assume a dominance of electrolytic conduction and associate the concentration with salt content, the tracer experiment can be well monitored by geoelectric measurements (Nguyen et al., 2009;Doetsch et al., 2012;Wagner et al., 2013). The fluid resistivity ρ f is obtained by a linear transformation: with a conductivity of groundwater σ 0 ¼ 0:1 mS=cm and a conversion factor after Sulzbacher et al. (2012) based on in-situ data. The Archie equation relates the bulk electrical resistivity of a medium ρ to the fluid resistivity ρ f and saturation S depending on porosity ϕ (Archie, 1942): The tortuosity factor a and the saturation exponent n are commonly assumed to be a ¼ 1 and n ¼ 2, respectively. The cementation exponent m is assumed to equal 2 (standard for sandstones). We assume no fluid saturation in the top layer and full saturation in the aquifer with a porosity of ϕ ¼ 0:3. To simulate synthetic data, i.e., apparent resistivities, we apply a dipole-dipole array with 41 equally spaced electrodes. Note that, compared to the Darcy simulation, ERT modelling requires specific boundary conditions and mesh refinement . Therefore, we create a suitable ERT forward mesh and interpolate the bulk resistivity values. From the 1600 advection time frames, we select 10 to simulate the ERT data by applying the resistivity values and the measuring scheme to the ERT Manager and call the simulate interface. The minimalist example code can be extended by the following code snippet: to create ERT dataset of apparent resistivities, which are shown in Fig. 7. The main characteristics of the tracer movement are clearly reflected by the measurement.
There are several approaches for inversion of time-lapse ERT data, like sequential inversion (e.g., Coscia et al., 2011), difference inversion (e.g., Bechtold et al., 2012), or 4-D inversion (Karaoulis et al., 2013). These approaches are available through pyGIMLi frameworks and result in a spatio-temporal resistivity distribution. Depending on temporal and spatial resolution, it might be possible to find out a single hydraulic conductivity value by analyzing breakthrough curves or similar techniques. However, in this case we go one step further and combine all simulation steps together into a sequentially coupled forward operator that maps a distribution of hydraulic conductivities (Fig. 4) into a set of apparent resistivities (Fig. 7) for each timestep. As a result, we are able to invert for the hydraulic conductivity directly (Mboh et al., 2012).
The forward operator provides the necessary interface and can be used by the basic inversion framework to perform a fully coupled hydrogeophysical inversion from measured apparent resistivities to an image of hydraulic conductivity within the aquifer. As an efficient and optimized exact generation of the Jacobian matrix is not easily possible for an arbitrarily coupled system, a brute-force Jacobian calculation is performed by solving the forward problem for each unknown model cell. As this task can be easily parallelized, it is automatically done if multiple processors are detected. Furthermore, we limit the number of unknown parameters regarding the possible resolution of the inverse problem and to reduce the numerical effort. We apply the pyGIMLi region concept to introduce three regions, of which two are treated as fixed regions with known hydraulic conductivities on the top and at the bottom of the domain. The assumed aquifer is subdivided into 256 rectangular parameter regions and further into a fine triangular forward mesh above to ensure numerical accuracy (Fig. 8). The inversion starts with a homogeneous model and converged fitting the simulated ERT data within the assumed data error. The hydraulic conductivity distribution illustrated by Fig. 8 resembles the main characteristics of the synthetic model both structurally and concerning the mean values.
Note that the example serves as a proof of concept and is intended to Fig. 7. Simulated dipole-dipole data of six consecutive measurements to monitor the migration of a saline tracer illustrated in Fig. 6. demonstrate that complex workflows can be easily realized using the flexible modelling and inversion framework presented. Yet it should be noted that geophysics-based reconstruction of hydraulic conductivity is considerably challenged if the scenario becomes more realistic (e.g., 3D, unsaturated flow, temperature effects).

Petrophysical joint inversion
Joint inversion of different geophysical techniques helps to improve both resolution and interpretability of the resulting images. Different data sets can be directly coupled, if there is a link to an underlying target parameter. This is also possible if the relations are not known beforehand (Heincke et al., 2017). The following example assumes the relations as known and demonstrates the framework concept for a petrophysically coupled joint inversion of two geophysical methods in a few simple steps. We use the existing method managers for ERT (pg.physics.ERTManager) and travel-time tomography (TT) (pg.physics.Refraction), which both provide ready-to-use simulation and inversion capabilities.
Initially, we apply both method managers to create synthetic data sets. Fig. 9 shows the used mesh (mMesh) for a laboratory-scale circular model (e.g., a sandstone column) with constant porosity ϕ ¼ 0:3 and heterogeneous water saturation S, which is investigated using ERT and ultrasonic tomography.
To create geophysical parameter distributions, we apply common empirical petrophysical models, e.g., Archie's equation (eq. (11)) that provides the resistivity ρ as a function of saturation S. Sonic velocity v (or its reverse, the slowness s) as a function of porosity ϕ and saturation S is given by the time-average equation (Wyllie et al., 1956): We assume a matrix velocity v m ¼ 4000 m/s, a water velocity v w ¼ 1484 m/s and an air velocity v a ¼ 343 m/s. pyGIMLi contains a collection of mechanisms for forward-and inverse value mapping provided by Trans classes. They can be used for data and model scaling in the inversion process, e.g., by using logarithmic barriers (Kim and Kim, 2011) to restrict parameters between given ranges. Similarly, petrophysical relations can be easily defined and also nested with other transforms as range constraints. Assume the intrinsic parameter p depends on the model parameter m. A transformation class needs to provide a forward and a backward transformation as well as the derivative of the transform, which is determined by a Newton algorithm if not specified.

Trans:fwdðSÞ
Trans:invðrhoÞ Trans:derivðrhoÞ ∂p ∂m We apply these classes and create transformations for the mentioned petrophysical relations under the names ArchieTrans for ERT and WyllieTrans for TT, respectively. The resistivity res and velocity vel are generated for a given saturation by running the simple code: Fig. 10 shows the images for the resulting mapped ρ (res) and v (vel) values. Note that although we show the commonly used velocity, the slowness is the intrinsic linear parameter used inside.
To create synthetic data sets, we assume 16 equally-spaced sensors on the circumferential boundary of the mesh. For the ERT modelling we build a complete dipole-dipole array. For the ultrasonic tomography we simulate the travel time for every possible sensor pair. The modelling itself is performed by calling the simulate interface command of the individual method managers, which create the data sets ertData for the resistivity and ttData for the velocity parameter distribution. Additionally, we add Gaussian noise with a standard deviation of 1% to the synthetic data sets.
To avoid an inverse crime, where identical meshes are used for both forward modelling and inversion (e.g. Colton and Kress, 1992), we create a new mesh for inversion, which holds no prior information on the anomalies to be reconstructed. Therefore, we create a second mesh pMesh representing the circular model domain that is coarser compared to the simulation mesh and contains no prior information on the anomalies. The simulated data are inverted on the new parameterization using the default inversion framework accessed through the interface command invert provided by both method managers and can be read as: Fig. 11 shows the resulting resistivity and velocity models. Generally, they recover the main structures of the synthetic model. The image of the resistivity distribution is more blurry but shows both the low and high resistive parts. In contrary, the velocity image shows sharper contrasts but lacks the low-velocity regions.
To reconstruct saturation values from the resulting images one can easily map the results back using the two transformation classes C. Rücker et al. Computers and Geosciences 109 (2017) 106-123 discussed above. However, we rather want a petrophysically constrained inversion as a more sophisticated way to directly recover the saturation b m ¼ fSg. To calculate the forward response F ðmÞ and the Jacobian J for a given saturation, we create a new forward operator called Petro-Modelling, which is method independent but initialized with preconfigured instances of a forward operator and a transformation class. It handles the necessary transformation steps and delegates the response and createJacobian commands appropriately. The forward response for the transformed problem is written as: The Jacobian matrix for the transformed problem is obtained by multiplication with the inner derivative of the transformation: For this we use a special matrix MultRightMatrix that holds an inner matrix (for which the Jacobian generation already exists) and a vector to be multiplied with from the right. The complete and concise implementation of the PetroModelling class is given in Appendix A and can be directly used with the generalized inversion framework.
As this is often used, an inversion framework PetroInversion for this type of petrophysical inversion is provided. It can be directly initialized by a method-specific manager and a transformation class. All relevant information is requested by the method managers provided and allows to construct a dedicated petrophysical forward operator with a corresponding inversion interface command (invert) and an additional argument limits for the expected model value ranges (log-transformation). The complete petrophysical inversion can then be written in two lines of code: Fig. 12 shows the resulting models for the petrophysically constrained inversion for saturation of both synthetic data sets.
The characteristics from both results is similar to the conventional  To combine the inherent information of both data sets in the estimation of the same parameter, a general framework for joint inversions is provided in pyGIMLi. The core for this framework is a forward operator called JointModelling, which allows a simple stacking of F different forward operators to create the appropriate model response and entries for the Jacobian matrix.
The model response is calculated by simply concatenating the responses of the individual operators. The Jacobian matrix for this type of concatenated problem can be efficiently formulated using a a special kind of matrix that we call BlockMatrix. This class allows arbitrary combinations of different matrix types that are directly used in the inverse solver. In this case (different data and different models) it has the shape of a block-Jacobi matrix with the individual Jacobian matrices on the diagonal. Block matrices is extremely helpful for a number of problems where independent models are simultaneously inverted, e.g., for LCI=SCI, time-lapse or spectral inversion . A minimal but complete implementation for the JointModelling class if given in Appendix B.
To keep the usage simple, we created a petrophysical joint inversion framework JointPetroInversion (see Appendix C), which combines the two forward operators PetroModelling and JointModelling. It is initialized by a set of method managers and corresponding petrophysical transforms. The framework manages the creation and delegation of the JointModelling forward operator, synchronizes data and error handling and provides the usual interface commands to perform the calculation via: Fig. 13 shows the resulting saturation image for the petrophysical joint inversion. The result of the joint inversion provides an improved image compared to both single inversion results as combines the benefits of both methods: It shows sharper contrasts than pure ERT and the regions of low and high saturation are more clearly visible compared to single TT.
In real life, the parameters of the petrophysical relations might not be known beforehand. However, it is straightforward to include a petrophysical parameter or even a distribution of it into the inversion process along with appropriate range constraints.

Conclusions
We have presented pyGIMLi, a versatile open-source framework for modelling and inversion in geophysics, which, due to its generalized and object-oriented design, is particularly useful to couple different measurement methods in joint or coupled inversions. This was demonstrated by a fully coupled inversion of time-lapse ERT data that allows to directly estimate the hydraulic conductivity distribution of the aquifer. The multiphysical forward operator provides the hydraulic and geoelectrical response of the tracer migration. It would thus be straightforward to also include hydrological data in the inversion to better constrain hydraulic aquifer properties.
A new petrophysical joint inversion approach was introduced and applied to ultrasonic and ERT combines the benefits of the two methods for a more accurate quantification of water saturation. In case of nonexactly known petrophysical parameters (e.g., porosity) this scheme could be easily extended, even without fixed relations. The availability of different forward operators and a generalized inversion framework  provides unprecedented means of mutual model constraints. For example, one could easily combine the joint inversion approach with the hydrogeophysical inversion so that a joint hydrogeophysical inversion for ERT data with hydrogeological (e.g. salinity or hydraulic head) data is achieved. All examples were generated with pyGIMLi version 1.0 and are fully reproducible using the scripts provided at http://cg17.pygimli.org.
Although the examples are based on simple 2D geometries for the sake of demonstration, we point out that the dimension of the problem is solely governed by the mesh provided, meaning that all code examples are directly transferable to more complex 3D geometries. The modular functionality is accessible from different abstraction levels and thereby offers entry points for a wide range of users with different degrees of programming experience: Application level Data owners can easily analyze, invert and visualize their data using predefined method managers.
Modelling level Users with custom forward operators can readily setup a corresponding inversion workflow with flexible discretization and regularization controls.

Equation level
Technically versatile practitioners can directly access the finite element and finite volume solvers to approach various physical problems.
There is a lot of potential for improving the package on all levels, starting from other PDE or element types over the integration of new methods to the extension of the inversion and regularization approaches.
Particularly, we see benefits of combination with other open Python/ Cþþ packages focusing on different methods.
To facilitate adoption and contribution of the geoscientific community, we have made the software package and all examples shown in this paper freely available under the permissive Apache 2.0 license. Corresponding downloads and a comprehensive documentation can be found on the project website http://www.pygimli.org. The source code is hosted on GitHub (https://github.com/gimli-org/gimli). By following modern software development principles such as unit testing and continuous integration, robustness and validity of existing and new functionality is ensured. It is anticipated that the presented software and ensuing developments will contribute to meet the modern data integration needs of researchers, practitioners, and students, particularly in, but not limited to, the hydrogeophysical community.

Acknowledgments
The constructive reviews of Rowan Cockett and two anonymous reviewers have considerably improved the manuscript. We thank Jennifer Geacone-Cruz, Nico Skibbe, and Thomas Heinze for careful reading and valuable hints. We would also like to express our appreciation for the growing number of pyGIMLi users and their important feedback enabling a continuous improvement of the software.

Appendices.
The following appendices contain the codes to carry out a petrophysical joint inversion as shown in 3.3. The two modelling operators for petrophysical modelling and joint modelling are used in the framework to carry out the inversion. Note that the implementations are independent on the actual methods and relations to be used.