1 
\chapter{Gravity Inversion}\label{Chp:cook:gravity inversion} 
2 

3 
\section{Introduction} 
4 

5 
In this part of the documentation we give an introduction on how to use the 
6 
\downunder module and the inversion driver functions to perform the inversion 
7 
of gravity and magnetic data. 
8 
The driver functions enable geologists and geophysicists to apply the 
9 
\downunder module quickly and in an easy way without requiring detailed 
10 
knowledge of the theory behind inversion or programming skills. 
11 
However, users who are interested in specializing or extending the inversion 
12 
capacity are referred to Part~\ref{part2} of this manual. 
13 
It is beyond the intention of this manual to give a detailed introduction to 
14 
geophysical inversion, in particular to the appropriate preprocessing of data 
15 
sets. 
16 

17 
The \downunder module described here is designed to calculate estimations for 
18 
the 3D distribution of density and/or susceptibility from 2D gravity and 
19 
magnetic data measured in ground or airborne surveys. 
20 
This process is generally called inversion of geophysical data. 
21 
Following the standard assumption it is assumed that the data are measured as 
22 
perturbation of an expected gravity and/or magnetic field of the Earth. 
23 
In this context measured gravity and magnetic data are in fact describing 
24 
anomalies in the gravity and magnetic field. 
25 
As a consequence the inversion process provides corrections to an average 
26 
density (typically $2670 kg/m^3$) and susceptibility (typically $0$). 
27 
So in the following we will always assume that given data are anomalies and 
28 
therefore not in all cases explicitly use the terms gravity anomalies or 
29 
density corrections but just use the terms gravity and density. 
30 

31 
In this chapter we will give a detailed introduction into usage of the driver 
32 
functions for inversion of gravity data. 
33 
In the following chapters~\ref{Chp:cook:magnetic inversion} and~\ref{Chp:cook:joint inversion} 
34 
we will discuss the inversion of magnetic data, and the joint inversion of 
35 
gravity and magnetic data using \downunder. 
36 
Note, that the principles introduced in this chapter apply to magnetic and 
37 
joint inversion so the presentation for these problem classes is kept short 
38 
and users interested in magnetic data only should still work through this 
39 
chapter on gravity data. 
40 

41 
To run the examples discussed you need to have \escript (version 3.3.1 or 
42 
newer) installed on your computer. 
43 
Moreover, if you want to visualize the results you need to have access to a 
44 
data plotting software which is able to process \VTK input files, e.g. 
45 
\mayavi or \VisIt. 
46 
As \mayavi can be easily obtained and installed for most platforms the 
47 
tutorial includes commands to visualize output files using \mayavi . 
48 
However, it is pointed out that \VisIt is the preferred visualization tool for 
49 
\escript as it can deal with very large data sets more efficiently. 
50 

51 
\begin{figure} 
52 
\centering 
53 
\includegraphics[width=0.7\textwidth]{QLDWestGravityDataPlot.png} 
54 
\caption{Gravity Anomaly Data in $mgal$ from Western Queensland, Australia 
55 
(see file \examplefile{data/QLDWest_grav.nc}). Data obtained from Geoscience Australia.} 
56 
%\AZADEH{Tracey R, Bacchin M, \& Wynne P. 2008. In preparation. AAGD07: A new absolute gravity datum for Australian gravity and new standards for the Australian National Gravity Database. Exploration Geophysics.}} 
57 
\label{FIG:P1:GRAV:0} 
58 
\end{figure} 
59 

60 
\begin{figure} 
61 
\centering% 
62 
\includegraphics[width=0.9\textwidth]{QLDGravContourMu10.png} 
63 
\caption{3D contour plot of the density distribution obtained by inversion of 
64 
file \file{data/QLDWest_grav.nc} (with $\mu=10$). 
65 
Colours represent values of density where high values are represented by 
66 
blue and low values are represented by red.} 
67 
\label{FIG:P1:GRAV:1} 
68 
\end{figure} 
69 

70 
\section{How does it work?} 
71 
The execution of the inversion is controlled by a script which, in essence, is 
72 
a text file and can be edited using any text editor. 
73 
The script contains a series of statements which are executed by an 
74 
interpreter which is an executable program reading the text file and executing 
75 
the statements linebyline. 
76 
In the case of \downunder the interpreter is \Python. 
77 
In order to be able to process the statements in each line of the script 
78 
certain rules (called syntax) need to be obeyed. 
79 
There is a large number of online tutorials for \Python available\footnote{e.g. 
80 
\url{http://www.tutorialspoint.com/python} and \url{http://doc.pyschools.com}}. 
81 
We also refer to the \escript cook book \cite{ESCRIPTCOOKBOOK} and user's 
82 
guide \cite{ESCRIPT} which is in particular useful for users who like to dive 
83 
deeper into \downunder. 
84 
For this part of the manual no \Python knowledge is required but it is 
85 
recommended that users acquire some basic knowledge on \Python as they 
86 
progress in their work with \downunder. 
87 

88 
The following script~\ref{code: gravity1}\footnote{The script is similar to 
89 
\examplefile{grav_netcdf.py} within the \escript example file directory.} is a 
90 
simple example to run an inversion for gravity data: 
91 

92 
\begin{pyc}\label{code: gravity1} 
93 
\ 
94 
\begin{python} 
95 
# Header: 
96 
from esys.downunder import * 
97 
from esys.weipa import * 
98 
from esys.escript import unitsSI as U 
99 

100 
# Step 1: set up domain 
101 
dom=DomainBuilder() 
102 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
103 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
104 
dom.fixDensityBelow(depth=40.*U.km) 
105 

106 
# Step 2: read gravity data 
107 
source0=NetCdfData(NetCdfData.GRAVITY, 'QLDWest_grav.nc', altitude=0.) 
108 
dom.addSource(source0) 
109 

110 
# Step 3: set up inversion 
111 
inv=GravityInversion() 
112 
inv.setSolverTolerance(1e4) 
113 
inv.setSolverMaxIterations(50) 
114 
inv.setup(dom) 
115 

116 
# Step 4: run inversion 
117 
inv.getCostFunction().setTradeOffFactorsModels(10.) 
118 
rho = inv.run() 
119 

120 
# Step 5: write reconstructed density to file 
121 
saveVTK("result.vtu", density=rho) 
122 
\end{python} 
123 
\end{pyc} 
124 
The result, in this case the density distribution, is written to an external 
125 
file for further processing. You can copy and paste the text of the script 
126 
into a file of any name, let's say for further reference we use the file 
127 
name \file{grav.py}. 
128 
It is recommended to use the extension \file{.py} to identify the file as a 
129 
\Python script. 
130 
We will discuss the statements of the script later in this chapter. 
131 

132 
Obviously the inversion needs to be fed with some gravity data. In this case 
133 
the data are loaded from the file \examplefile{QLDWest_grav.nc} which is 
134 
available in the \escript example file directory. 
135 
The data in this example are given in the \netcdf file format. 
136 
After you have copied this file into the directory in which you have saved the 
137 
script \file{grav.py} you can run the program using the command line 
138 
\begin{verbatim} 
139 
runescript grav.py 
140 
\end{verbatim} 
141 
We are running \file{grav.py} through the \escript startup command since 
142 
\escript is used as a back end for the inversion algorithm\footnote{Please see 
143 
the \escript user's guide~\cite{ESCRIPT} on how to run your script in parallel 
144 
using threading and/or MPI.}. 
145 
Obviously it is assumed that you have an installation of \escript available on 
146 
your computer, see \url{https://launchpad.net/escriptfinley}. 
147 

148 
After the execution has successfully completed you will find the result file 
149 
\file{result.vtu} in the directory from where you have started the execution 
150 
of the script. 
151 
The file has the \VTK format and can be imported easily into many 
152 
visualization tools. 
153 
One option is the \mayavi package which is available on most platforms. 
154 
You can invoke the the visualization using the commands 
155 
\begin{verbatim} 
156 
mayavi2 d result.vtu m SurfaceMap 
157 
\end{verbatim} 
158 
from the command line. 
159 
Figure~\ref{FIG:P1:GRAV:1} shows the result of this inversion as a contour 
160 
plot\footnote{These plots were generated by \VisIt.}, while the gravity 
161 
anomaly data is shown in Figure~\ref{FIG:P1:GRAV:0}. 
162 
We will discuss data later in Section~\ref{SEC:P1:GRAV:DATA}. 
163 

164 
Let us have a closer look at the script\footnote{In \Python lines starting 
165 
with `\#` are comments and are not processed further.}: Besides the header 
166 
section one can separate the script into five steps: 
167 
\begin{enumerate} 
168 
\item set up domain on which the inversion problem is solved 
169 
\item load the data 
170 
\item setup the inversion problem 
171 
\item run the inversion 
172 
\item further processing of the result. Here we write the reconstructed 
173 
density distribution to a file. 
174 
\end{enumerate} 
175 
In the following we will discuss the steps of the scripts in more detail. 
176 
Before we do this it is pointed out that the header section, following 
177 
\Python conventions, makes all required packages available to access within 
178 
the script. 
179 
At this point we will not discuss this in more details but emphasize that the 
180 
header section is a vital part of the script. 
181 
It is is required in each \downunder inversion script and should not be 
182 
altered except if additional modules are needed. 
183 

184 
\begin{figure} 
185 
\centering 
186 
\includegraphics[width=\textwidth]{dom2D.pdf} 
187 
\caption{2D domain setup for gravity inversion} 
188 
\label{FIG:P1:GRAV:2} 
189 
\end{figure} 
190 

191 
\section{Creating the Inversion Domain} 
192 
First step in script~\ref{code: gravity1} is the definition of the domain over 
193 
which the inversion is performed. 
194 
We think of the domain as a block with orthogonal, plain faces. 
195 
Figure~\ref{FIG:P1:GRAV:2} shows the setup for a twodimensional domain 
196 
(see also Figure~\ref{fig:domainBuilder} for 3D). 
197 
The lateral coordinates along the surface of the Earth are denoted by $x$ and 
198 
$y$ (only $x$direction is used in 2D). 
199 
The $z$ direction defines the vertical direction where the part above the 
200 
surface has positive coordinates and the subsurface negative coordinates. 
201 
The height of the section above the surface, which is assumed to be filled 
202 
with air, needs to be set by the user. 
203 
The inversion assumes that the density in the section is known to be 
204 
zero\footnote{Always keeping in mind that these are not absolute values but 
205 
anomalies.}. 
206 
The density below the surface is unknown and is calculated through the 
207 
inversion. The user needs to specify the depth below the surface in which the 
208 
density is to be calculated. 
209 
The lateral extension of the domain is defined by the data sets fed into the 
210 
inversion. 
211 
It is chosen large enough to cover all data sets (in case more than one is 
212 
used). In order to reduce the impact of the boundary a padding zone around the 
213 
data sets can be introduced. 
214 

215 
\begin{figure} 
216 
\centering 
217 
\includegraphics[width=\textwidth]{dom2DFEM.pdf} 
218 
\caption{Cell distribution and boundary conditions for a 2D domain} 
219 
\label{FIG:P1:GRAV:3} 
220 
\end{figure} 
221 

222 
The reconstruction of the gravity field from an estimated density distribution 
223 
is the key component of the inversion process. 
224 
To do this \downunder uses the finite element method (FEM). 
225 
We need to introduce a grid over the domain, see Figure~\ref{FIG:P1:GRAV:3}. 
226 
The number of vertical cells is set by the user while the number of horizontal 
227 
cells is derived from the grid spacing of the gravity data set(s). 
228 
It is assumed that gravity field data given are constant across a cell. 
229 
To be able to reconstruct the gravity field some assumptions on the values of 
230 
the gravity field on the domain boundary have to be made. 
231 
\downunder assumes that on all faces the lateral gravity field component 
232 
equals zero. No assumptions on the horizontal components are 
233 
made\footnote{It is assumed that the gravity potential equals zero on the top 
234 
and bottom surface, see Section~\ref{sec:forward gravity} for details}% 
235 
\footnote{Most inversion codes use Green's functions over an unbounded domain 
236 
to reconstruct the gravity field. This approach makes the assumption that the 
237 
gravity field (or potential) is converging to zero when moving away from the 
238 
region of interest. The boundary conditions used here are stronger in the 
239 
sense that the lateral gravity component is enforced to be zero in a defined 
240 
distance of the region of interest but weaker in the sense that no constraint 
241 
on the horizontal component is applied.}. 
242 

243 
In script~\ref{code: gravity1} the statement 
244 
\begin{verbatim} 
245 
dom=DomainBuilder() 
246 
\end{verbatim} 
247 
creates something like a factory to build a domain. 
248 
We then define the features of the domain we would like to create: 
249 
\begin{verbatim} 
250 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
251 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
252 
\end{verbatim} 
253 
Here we specify the depth of the domain to $40 km$, the thickness of the air 
254 
layer above the surface to $6 km$ and the number of vertical cells to $25$. 
255 
We also introduce a lateral padding of $20 \%$ of the expansion of the gravity 
256 
data on each side of the data and in both lateral directions. 
257 

258 
In some cases it can be appropriate to assume that the density below a certain 
259 
depth is zero\footnote{As we are in fact calculating density corrections this 
260 
means that the density is assumed to be the average density.}. 
261 
The statement 
262 
\begin{verbatim} 
263 
dom.fixDensityBelow(depth=40.*U.km) 
264 
\end{verbatim} 
265 
introduces this constraint. 
266 
As in the case discussed here if the depth for zero density is not less than 
267 
the depth of the domain no constraint at depth is applied to the density. 
268 

269 
\downunder uses the metrekilogramsecond based International System of Units 
270 
(SI)~\footnote{see \url{http://en.wikipedia.org/wiki/International_System_of_Units}}\index{SI}. 
271 
So all values must be converted to appropriate units. 
272 
This does not apply to geographic coordinates which in \downunder are given in 
273 
fractional degrees (as a floating point number) to represent longitude and 
274 
latitude. In the script we have used the expression 
275 
\begin{verbatim} 
276 
depth=40.*U.km 
277 
\end{verbatim} 
278 
to define the depth of the domain to $40 km$. 
279 
The expression \verbU.km denotes the unit $km$ (kilometer) and ensures 
280 
appropriate conversion of the value $40$ into the base unit $m$ (meter). 
281 
It is recommended to add units to values (where present) in order to make sure 
282 
that the final values handed into \downunder is given with the appropriate 
283 
units. 
284 
The physical units module of \escript, which we have imported here under the 
285 
name \verbU in the script header, defines a large number of physical units 
286 
and constants, please see~\cite{ESCRIPT} and~\cite{ESCRIPTONLINE}. 
287 

288 
\section{Loading Gravitational Data}\label{SEC:P1:GRAV:DATA} 
289 
In practice gravity acceleration is measured in various ways, for instance by 
290 
airborne surveys~\cite{Telford1990a}. 
291 
\downunder assumes that all data supplied as input are already appropriately 
292 
preprocessed. In particular, corrections for 
293 
\begin{itemize} 
294 
\item freeair, to remove effects from altitude above ground; 
295 
\item latitude, to remove effects from ellipsoidicity of the Earth; 
296 
\item terrain, to remove effects from topography 
297 
\end{itemize} 
298 
must have been applied to the data. 
299 
In general, data prepared in such a form are called Bouguer anomalies~\cite{Telford1990a}. 
300 

301 
To load gravity data into \downunder the data are given on a plane parallel 
302 
to the surface of the Earth at a constant altitude, see 
303 
diagram~\ref{FIG:P1:GRAV:2}. 
304 
The data need to be defined over a rectangular grid covering a subregion of 
305 
the Earth surface. 
306 
The grid uses a geographic coordinate system with latitudes and longitudes 
307 
assumed to be given in the Clarke 1866 geodetic system. 
308 
Figure~\ref{FIG:P1:GRAV:0} shows an example of such a data set from the 
309 
western parts of Queensland, Australia. 
310 
The data set covers a rectangular region between $140^o$ and $141^o$ east 
311 
and between $20^o$ and $21^o$ south. 
312 
Notice that latitude varies between $90^o$ to $90^o$ where negative signs 
313 
refer to places in the southern hemisphere and longitude varies between 
314 
$180^o$ to $180^o$ where negative signs refer to places in the western 
315 
hemisphere. 
316 
The colour at a location represents the value of the vertical Bouguer gravity 
317 
anomaly at this point at the surface of the Earth. 
318 
Values in this data set range from $160 \; mgal$ to about $500 \; mgal$\footnote{The unit 
319 
$mgal$ means milli $gal$ (galileo) with $1 \; gal = 0.01 \frac{m}{sec^2}$.} 
320 
over a $121 \times 121$ grid. 
321 

322 
In general, a patch of gravity data need to be defined over a plane 
323 
\verbNX $\times$ \verbNY where \verbNX and \verbNY define the number 
324 
of grid lines in the longitude (\verbX) and the latitude (\verbY) 
325 
direction, respectively. 
326 
The grid is spanned from an origin with spacing \verbDELTA_X and 
327 
\verbDELTA_Y in the longitude and the latitude direction, respectively. 
328 
The gravity data for all grid points need to be given as an \verbNX 
329 
$\times$ \verbNY array. 
330 
If available, measurement errors can be associated with the gravity data. 
331 
The values are given as an \verbNX $\times$ \verbNY array matching the 
332 
shape of the gravity array. 
333 
Note that data need not be available on every single point of the grid, see 
334 
Section~\ref{SEC:P1:GRAV:REMARK:DATAHOLES} for more information on this. 
335 

336 
Currently, two data file formats are supported, namely \emph{ER Mapper Raster} 
337 
\cite{ERMAPPER} files and \netcdf~\cite{NETCDF} files. 
338 
In the examples of this chapter we use \netcdf files and refer to 
339 
Section~\ref{SEC:P1:GRAV:REMARK:ERMAPPER} and 
340 
Section~\ref{sec:ref:DataSource:ERM} for more information on using ER Mapper 
341 
Raster files. 
342 
If you have data in any other format you have the option of writing a suitable 
343 
reader (for advanced users, see Chapter~\ref{Chp:ref:data sources}) or, 
344 
assuming you are able to read the data in \Python, refer to the example 
345 
script \examplefile{create_netcdf.py} which shows how to create a file 
346 
in the \netcdf file format~\cite{NETCDF} compatible with \downunder from 
347 
a data array. 
348 

349 
\CIHAN{Use coarse data 10x10 or so?} 
350 
In script~\ref{code: gravity1} we use the statement 
351 
\begin{verbatim} 
352 
source0=NetCdfData(NetCdfData.GRAVITY, 'QLDWest_grav.nc', altitude=0.) 
353 
\end{verbatim} 
354 
to load the gravity data stored in \examplefile{QLDWest_grav.nc} in the \netcdf format. 
355 
Within the script the data set is now available under the name \verbsource0. 
356 
We need to link the data set to the \verbDomainBuilder using 
357 
\begin{verbatim} 
358 
dom.addSource(source0) 
359 
\end{verbatim} 
360 
At the time the domain for the inversion is built the \verbDomainBuilder 
361 
will use the information about origin, extent and spacing along with the other 
362 
options provided to build an appropriate domain. 
363 
As at this point a flat Earth is assumed geographic coordinates used to 
364 
represent data in the input file are mapped to a (local) Cartesian coordinate 
365 
system. This is achieved by projecting the geographic coordinates into the 
366 
\emph{Universal Transverse Mercator} (UTM) coordinate system\footnote{See e.g. 
367 
\url{http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system}}. 
368 

369 
It is possible to add several, possibly overlapping data sets to the same 
370 
domain: 
371 
\begin{verbatim} 
372 
source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc', altitude=0.) 
373 
source1=NetCdfData(NetCdfData.GRAVITY, 'data1.nc', altitude=0.) 
374 
dom.addSource(source0) 
375 
dom.addSource(source1) 
376 
\end{verbatim} 
377 
However, due to the coordinate transformation all data sets must be located in 
378 
the same UTM zone. If a single dataset crosses UTM zones only the zone of the 
379 
central longitude is used when projecting. 
380 
For example, if one data set lies mostly in zone 51 but contains areas of zone 
381 
52, it is transformed using zone 51. 
382 
In this case more data from zone 51 can be added, but not from any other zone. 
383 

384 
It is important to keep an eye on the complexity of the inversion. 
385 
A good measure is the total number of cells being used. 
386 
Assume we have given a data set on a $20 \times 20$ grid and we add lateral 
387 
padding of, say, $20 \%$ to each side of the data, the lateral number of cells 
388 
becomes $(20\cdot 1.4)\times (20\cdot 1.4)=1.4^2\cdot 20^2\approx 2\cdot 10^2=800$. 
389 
If we use $20$ cells in the vertical direction we end up with a total number 
390 
of $800 \times 20 = 16,000$ cells. 
391 
This size can be easily handled by a modern desktop PC. 
392 
If we increase the grid size of the data to $40 \times 40$ points and use $40$ 
393 
cells in the vertical extent we get a total of $(2\cdot 40^2)\cdot 40=128,000$ 
394 
cells, a problem size which is considerably larger but can still be handled by 
395 
a desktop computer. 
396 
Taking this one step further, if the amount of data is increased to 
397 
$200\times 200$ points and we use $200$ cells in the vertical extent the 
398 
domain will contain $16,000,000$ ($16$ million) cells. 
399 
This scenario requires a computer with enough memory and (a) fast processor(s) 
400 
to run the inversion. 
401 
This estimate of complexity growth applies to the case where the increase 
402 
of data grid size is driven by an increase of resolution where it is 
403 
recommended to increase the vertical resolution in synch with the lateral 
404 
resolution. Note that if more than one data set is used the target resolution 
405 
will be the resolution of the finest data set. 
406 
In other cases the expansion of the region of interest drives an increase of 
407 
data grid size and the increase of total number of cells is less dramatic as 
408 
the vertical number of cells can remain constant while keeping a balanced 
409 
resolution in vertical and lateral direction. 
410 

411 
\section{Setting up the Inversion and Running it} 
412 
We are now at step three of script~\ref{code: gravity1} in which the actual 
413 
inversion is set up. 
414 
First we create an empty inversion under the name \verbinv: 
415 
\begin{verbatim} 
416 
inv=GravityInversion() 
417 
\end{verbatim} 
418 
As indicated by the name we can use \verbinv to perform an inversion of 
419 
gravity data\footnote{\verbGravityInversion is a driver with a simplified 
420 
interface which is provided for convenience. See Part~\ref{part2} for more 
421 
details on how to write inversion scripts with more general functionality, e.g. 
422 
constraints.}. The inversion is an iterative process which sequentially 
423 
calculates updates to the density distribution in an attempt to improve the 
424 
match of the gravity field produced by the density distribution with the data. 
425 
Termination of the iteration is controlled by the tolerance which is set by 
426 
the user: 
427 
\begin{verbatim} 
428 
inv.setSolverTolerance(1e4) 
429 
\end{verbatim} 
430 
Here we set the tolerance to $10^{4}$, i.e. the iteration is terminated if 
431 
the maximum density correction is less than or equal to $10^{4}$ relative to 
432 
the maximum value of estimated density anomaly. 
433 
In case the iteration does not converge a maximum number of iteration steps is 
434 
set: 
435 
\begin{verbatim} 
436 
inv.setSolverMaxIterations(50) 
437 
\end{verbatim} 
438 
If the maximum number of iteration steps (here $50$) is reached the iteration 
439 
process is aborted and an error message is printed. 
440 
In this case you can try to rerun the inversion with a larger value for the 
441 
maximum number of iteration steps. 
442 
If even for a very large number of iteration steps no convergence is achieved, 
443 
this is a strong indicator that the inversion has not been set up properly. 
444 

445 
The statement 
446 
\begin{verbatim} 
447 
inv.setup(dom) 
448 
\end{verbatim} 
449 
links the inversion with the domain and the data. 
450 
At this step  as we are solving a gravity inversion problem  only 
451 
gravitational data attached to the domain builder \verbdom are considered. 
452 
Internally a cost function $J$ is created which is minimized during the 
453 
inversion iteration. 
454 
It is a combination of a measure of the data misfit of the gravity field from 
455 
the given density distribution and a measure of the smoothness of the density 
456 
distribution. 
457 
The latter is often called the regularization term. 
458 
By default the gradient of density is used as the regularization term, see 
459 
also Section~\ref{SEC:P1:GRAV:REMARK:REG}. 
460 
Obviously, the result of the inversion is sensitive to the weighting between 
461 
the misfit and the regularization. 
462 
This tradeoff factor $\mu$ for the misfit function is set by the following 
463 
statement: 
464 
\begin{verbatim} 
465 
inv.getCostFunction().setTradeOffFactorsModels(0.1) 
466 
\end{verbatim} 
467 
Here we set $\mu=0.1$. The statement \verbinv.setup must appear in the 
468 
script before setting the tradeoff factor. 
469 
A small value for the tradeoff factor $\mu$ will give more emphasis to the 
470 
regularization component and create a smoother density distribution. 
471 
A large value of the tradeoff factor $\mu$ will emphasize the misfit more 
472 
and typically creates a better fit to the data and a rougher density 
473 
distribution. 
474 
It is important to keep in mind that the regularization reduces noise in the 
475 
date and in fact gives the problem a unique solution. 
476 
Consequently, the tradeoff factor $\mu$ may not be chosen too large in order 
477 
control the noise on the solution and ensure convergence in the iteration 
478 
process. 
479 

480 
We can now actually run the inversion: 
481 
\begin{verbatim} 
482 
rho = inv.run() 
483 
\end{verbatim} 
484 
The answer as calculated during the inversion is returned and can be accessed 
485 
under the name \verbrho. 
486 
As pointed out earlier the iteration process may fail in which case the 
487 
execution of the script is aborted with an error message. 
488 

489 
\section{Taking a Look} 
490 
In the final step of script~\ref{code: gravity1} the calculated density 
491 
distribution is written to an external file. 
492 
A popular file format used by several visualization packages such as 
493 
\VisIt~\cite{VISIT} and \mayavi~\cite{MAYAVI} is the \VTK file format. 
494 
The result of the inversion which has been named \verbrho can be written to 
495 
the file \file{result.vtu} by adding the statement 
496 
\begin{verbatim} 
497 
saveVTK("result.vtu", density=rho) 
498 
\end{verbatim} 
499 
at the end of script. 
500 
The inversion solution is tagged with the name \verbdensity in the result 
501 
file, however any other name for the tag could be used. 
502 
As the format is textbased (as opposed to binary) \VTK files tend to be very 
503 
large and take compute time to create, in particular when it comes to large 
504 
numbers of cells ($>10^6$). 
505 
For large problems it is more efficient to use the \SILO file format~\cite{SILO}. 
506 
\SILO files tend to be smaller and are faster generated and read. 
507 
It is the preferred format to import results into the visualization program 
508 
\VisIt~\cite{VISIT} which is particularly suited for the visualization of 
509 
large data sets. 
510 
Inversion results can directly exported into \SILO files using the statement 
511 
\begin{verbatim} 
512 
saveSilo("result.silo", density=rho) 
513 
\end{verbatim} 
514 
replacing the \verbsaveVTK(...) statement. 
515 
Similar to \VTK files the result \verbrho is tagged with the name 
516 
\verbdensity so it can be identified in the visualization program. 
517 

518 
\begin{figure} 
519 
\begin{center} 
520 
\subfigure[$\mu=0.1$]{% 
521 
\label{FIG:P1:GRAV:10 MU01} 
522 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu01.png} 
523 
}% 
524 
\subfigure[$\mu=1.$]{% 
525 
\label{FIG:P1:GRAV:10 MU1} 
526 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu1.png} 
527 
}\\ %  End of the first row % 
528 
\subfigure[$\mu=10.$]{% 
529 
\label{FIG:P1:GRAV:10 MU10} 
530 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu10.png} 
531 
}% 
532 
\subfigure[$\mu=100.$]{% 
533 
\label{FIG:P1:GRAV:10 MU100} 
534 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu100.png} 
535 
}\\ %  End of the second row % 
536 
\subfigure[$\mu=1000.$]{% 
537 
\label{FIG:P1:GRAV:10 MU1000} 
538 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu1000.png} 
539 
}% 
540 
\end{center} 
541 
\caption{3D contour plots of gravity inversion results with data from 
542 
Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
543 
factor $\mu$. Visualization has been performed in \VisIt. 
544 
\AZADEH{check images.}} 
545 
\label{FIG:P1:GRAV:10} 
546 
\end{figure} 
547 

548 
\begin{figure} 
549 
\begin{center} 
550 
\subfigure[$\mu=0.1$]{% 
551 
\label{FIG:P1:GRAV:11 MU01} 
552 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu01.png} 
553 
}% 
554 
\subfigure[$\mu=1.$]{% 
555 
\label{FIG:P1:GRAV:11 MU1} 
556 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu1.png} 
557 
}\\ %  End of the first row % 
558 
\subfigure[$\mu=10.$]{% 
559 
\label{FIG:P1:GRAV:11 MU10} 
560 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu10.png} 
561 
}% 
562 
\subfigure[$\mu=100.$]{% 
563 
\label{FIG:P1:GRAV:11 MU100} 
564 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu100.png} 
565 
}\\ %  End of the second row % 
566 
\subfigure[$\mu=1000.$]{% 
567 
\label{FIG:P1:GRAV:11 MU1000} 
568 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu1000.png} 
569 
}% 
570 
\end{center} 
571 
\caption{3D slice plots of gravity inversion results with data from 
572 
Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
573 
factor $\mu$. Visualization has been performed \VisIt. 
574 
\AZADEH{check images.}} 
575 
\label{FIG:P1:GRAV:11} 
576 
\end{figure} 
577 

578 
Figures~\ref{FIG:P1:GRAV:10} and~\ref{FIG:P1:GRAV:11} show two different 
579 
styles of visualization generated in \VisIt using the result of the inversion 
580 
of the gravity anomalies shown in Figure~\ref{FIG:P1:GRAV:0}. 
581 
The inversions have been performed with different values for the model 
582 
tradeoff factor $\mu$. 
583 
The visualization shows clearly the smoothing effect of lower values for the 
584 
tradeoff factors. 
585 
For larger values of the tradeoff factor the density distribution becomes 
586 
rougher showing larger details. 
587 
Computational costs are significantly higher for larger tradeoff factors. 
588 
Moreover, noise in the data has an higher impact on the result. 
589 
Typically several runs are required to adjust the value for the tradeoff 
590 
factor to the datasets used. 
591 

592 
For some analysis tools it is useful to process the results in form of Comma 
593 
Separated Values (\CSV)~\footnote{see 
594 
\url{http://en.wikipedia.org/wiki/Commaseparated_values}.}. 
595 
Such a file can be created using the statement 
596 
\begin{verbatim} 
597 
saveDataCSV("result.csv", x=rho.getFunctionSpace().getX(), density=rho) 
598 
\end{verbatim} 
599 
in the script. 
600 
This will create a \file{result.csv} with columns separated by a comma. 
601 
Each row contains the value of the density distribution and the three 
602 
coordinates of the corresponding location in the domain. 
603 
There is a header specifying the meaning of the corresponding column. 
604 
\LG{Add example output CSV}. 
605 
Notice that rows are not written in a particular order and therefore, if 
606 
necessary, the user has to apply appropriate sorting of the rows. 
607 
Columns are written in alphabetic order of their corresponding tag names. 
608 
For the interested reader: the statement \verbrho.getFunctionSpace() returns 
609 
the type used to store the density data \verbrho. 
610 
The \verbgetX() method returns the coordinates of the sampling points used 
611 
for the particular type of representation, see~\cite{ESCRIPT} for details. 
612 

613 

614 
\section{Remarks} 
615 

616 
\subsection{Data With Holes}\label{SEC:P1:GRAV:REMARK:DATAHOLES} 
617 
\CIHAN{ADD some remarks holes in the data, example?} 
618 

619 

620 
\subsection{ER Mapper Raster File}\label{SEC:P1:GRAV:REMARK:ERMAPPER} 
621 
\CIHAN{ADD an example how to use ER Mapper raster files} 
622 

623 
\subsection{Regularization Term}\label{SEC:P1:GRAV:REMARK:REG} 
624 
The \verbGravityInversion class supports the following form for the 
625 
regularization: 
626 
\begin{equation} 
627 
\int w^{(0)} \cdot \rho^2 + w^{(1)}_0 \rho_{,0}^2 + w^{(1)}_1 \rho_{,1}^2 + w^{(1)}_2 \rho_{,2}^2\; dx 
628 
\end{equation} 
629 
where the integral is calculated across the entire domain. 
630 
$\rho$ represents the density distribution where $\rho_{,0}$ $\rho_{,1}$ and 
631 
$\rho_{,2}$ are the spatial derivatives of $\rho$ with respect to the two 
632 
lateral and the vertical direction, respectively. 
633 
$w^{(0)}$, $w^{(1)}_0$, $w^{(1)}_1$ and $w^{(1)}_2$ are weighting 
634 
factors\footnote{A more general form, e.g. spatially variable values for the 
635 
weighting factors, is supported, see Part~\ref{part2}}. 
636 
By default these are $w^{(0)}=0$, $w^{(1)}_0=w^{(1)}_1=w^{(1)}_2=1$. 
637 
Other weighting factors can be set in the inversion setup. 
638 
For instance to set $w^{(0)}=10$, $w^{(1)}_0=w^{(1)}_1=0$ and $w^{(1)}_2=100$ 
639 
use the statement: 
640 
\begin{verbatim} 
641 
inv.setup(dom, w0=10, w1=[0,0,100]) 
642 
\end{verbatim} 
643 
It is pointed out that the weighting factors are rescaled in order to improve 
644 
numerical stability. Therefore the relative size of the weighting factors is 
645 
relevant and using 
646 
\begin{verbatim} 
647 
inv.setup(dom, w0=0.1, w1=[0,0,1]) 
648 
\end{verbatim} 
649 
would lead to the same regularization as the statement above. 
650 

651 
% 
652 
% \section{Reference} 
653 
% 
654 
% There are three examples for 2D and 3D gravity inversions with artificial input data. 
655 
% In first step, an area with synthetic density section was suggested. Then based on forward modeling its gravitational data was collected. Afterwards with generated gravity data, escripts simulated a volume of inverted density. Whilst new density mass could be compared with the synthetic density section to verify the inversion. 
656 
% 
657 
% Some of the presumptions were the same for all of the examples to simplify the situation to make a logical comparison between synthetic input and output. which is as followed: 
658 
% 
659 
% \begin{verbatim} 
660 
% n_cells_in_data=100 
661 
% depth_offset=0.*U.km 
662 
% l_data = 100 * U.km 
663 
% l_pad=40*U.km 
664 
% THICKNESS=20.*U.km 
665 
% l_air=6*U.km 
666 
% \end{verbatim} 
667 
% 
668 
% The others assumptions comes with each examples. 
669 
% \begin{enumerate} 
670 
% \item A 2D density section with a maximum in center was assumed. The reference density and inverted will be shown. The padding area is excluded. (\ref{fig:gravity2D1}) 
671 
% 
672 
% \begin{verbatim} 
673 
% n_cells_in_data=100 
674 
% n_humbs_h= 1 
675 
% n_humbs_v=1 
676 
% mu=100 
677 
% \end{verbatim} 
678 
% 
679 
% \begin{figure} 
680 
% \centering 
681 
% \includegraphics[width=\textwidth]{grav2D1.png} 
682 
% \caption{2D density model up) reference down) result} 
683 
% \label{fig:gravity2D1} 
684 
% \end{figure} 
685 
% 
686 
% 
687 
% \item A 2D density properties with two maximum in corners and one minimum in the center was inverted. The result have eliminated the effects in padding area.(\ref{fig:gravity2D3}) 
688 
% 
689 
% \begin{verbatim} 
690 
% n_cells_in_data=100 
691 
% n_humbs_h= 3 
692 
% n_humbs_v=1 
693 
% mu=100 
694 
% \end{verbatim} 
695 
% 
696 
% \begin{figure} 
697 
% \centering 
698 
% \includegraphics[width=\textwidth]{grav2D3.png} 
699 
% \caption{2D density model up) reference down) result} 
700 
% \label{fig:gravity2D3} 
701 
% \end{figure} 
702 
% 
703 
% \item A 3D model with a maximum in the center was used as input data and the result after simulation in shown in the next image which determined a good distribution of density through the model in the main area.(\ref{fig:gravity3D} and \ref{fig:gravity3D1}) 
704 
% 
705 
% \begin{verbatim} 
706 
% n_cells_in_data=50 
707 
% n_humbs_h= 1 
708 
% n_humbs_v=1 
709 
% mu=10 
710 
% \end{verbatim} 
711 
% 
712 
% \begin{figure} 
713 
% \centering 
714 
% \includegraphics[width=\textwidth]{density3Dref.png} 
715 
% \caption{3D density model of reference as synthetic data} 
716 
% \label{fig:gravity3D} 
717 
% \end{figure} 
718 
% 
719 
% \begin{figure} 
720 
% \centering 
721 
% \includegraphics[width=\textwidth]{gravity3D.png} 
722 
% \caption{3D density model of result} 
723 
% \label{fig:gravity3D1} 
724 
% \end{figure} 
725 
% \end{enumerate} 
726 
