Computational methods in hydrogeophysicsbyKlara SteklovaB. Applied Ecology, Czech University of Life Science, Prague, 2006M. Ecohydrology and Groundwater management, Wageningen University, 2009a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoral studies(Geophysics)The University of British Columbia(Vancouver)March 2017c© Klara Steklova, 2017AbstractParameter and state estimation for groundwater models within a coupled hydrogeophysicalframework has become common in the last few years as it has been shown that suchestimates are usually better than those from a single data inversion. Different approacheshave been suggested in literature to combine the essentially two different modalities in orderto obtain better estimates for groundwater models, and improve monitoring of processessuch as solute transport. However, the coupled approaches usually come at a price ofhigher computational cost and difficulties in coupling the geophysical and groundwaterinverse problems.Unlike in other studies, we developed both the groundwater and geophysical modelsin the same computational environment in order to test different minimization strategies.When solving the coupled inverse problem, the objective function consists of data misfit andregularization terms as well as a coupling term that relates groundwater and geophysicalstates. We present a novel approach to solve the inverse problem using an AlternatingDirection Method of Multipliers (ADMM) to minimize the coupled objective function.ADMM enables us to treat the groundwater and geophysical part separately and thus useexisting software with minor changes.However, ADMM as well as many other coupled approaches relies on implementingsome petrophysical relationship to couple the groundwater and geophysical variable. Suchrelationships are usually uncertain and hard to parametrize for a large region and canpotentially produce solute mass errors in the final model estimates. Therefore, in this thesiswe examine coupled approaches that replace the fixed petrophysical relationship by a moreloose structure similarity constraint. Besides, we propose efficient computational methodsto minimize the objective function when there is no explicit petrophysical constraint.All approaches were tested on 3D synthetic examples. In the solute tracer test we es-timated hydraulic conductivity or solute distribution using a structure coupled inversion,iiand were able to reduce the errors compared to a single data inversion alone. For a morecomplex example of seawater intrusion we implemented the ADMM method, and obtainedbetter estimates for the solute distribution compared to just considering each data sepa-rately, or solving the problem with a simple coupled approach.iiiPrefaceThis dissertation is my original work made under supervision of Dr. Haber during thestudy program at the University of British Columbia.The first three chapters review the current modeling practices and state of the art inhydrogeophysics. Chapter 4 presents new ideas for the joint inversion minimization andstructure coupled inversion which came from conversations with Dr. Haber.The code implementation and numerical experiments were carried out by myself withDr. Haber’s support and advice, and form the Chapter 5. The results in the first sectionof chapter 5 were published in [91], and the results in the second section are part of themanuscript, which is currently under review [90].ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Geophysical data in hydrology . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Hydrogeophysics: different approaches to couple the modalities . . . . . . . 51.3 Coupled hydrogeophysical inversion: The state of art . . . . . . . . . . . . . 71.3.1 Structure - coupled inversion . . . . . . . . . . . . . . . . . . . . . . 91.4 Thesis objectives and research contributions . . . . . . . . . . . . . . . . . . 112 Forward modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1 Groundwater model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1.1 Modeling of solute tracer movement . . . . . . . . . . . . . . . . . . 152.1.2 Variable density flow (VDF) model . . . . . . . . . . . . . . . . . . . 212.1.3 Testing the GW models . . . . . . . . . . . . . . . . . . . . . . . . . 242.2 Geophysical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.2.1 Direct current resistivity survey . . . . . . . . . . . . . . . . . . . . . 272.2.2 Time domain electromagnetic survey (TDEM) . . . . . . . . . . . . 28v3 Inverse problem for a single model . . . . . . . . . . . . . . . . . . . . . . 323.1 Ill-posedness and regularization . . . . . . . . . . . . . . . . . . . . . . . . . 323.2 Gauss-Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3 Sensitivities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364 Coupled inverse problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.1 Coupling the GW and geophysical variables . . . . . . . . . . . . . . . . . . 424.1.1 Petrophysical relationships . . . . . . . . . . . . . . . . . . . . . . . 424.1.2 Structure-based coupling . . . . . . . . . . . . . . . . . . . . . . . . 444.1.3 Discretization of the coupling terms . . . . . . . . . . . . . . . . . . 464.2 Formulation of the coupled problem . . . . . . . . . . . . . . . . . . . . . . 484.3 Minimization with an explicit petrophysical constraint . . . . . . . . . . . . 494.3.1 ADMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.4 Minimization without an explicit petrophysical constraint . . . . . . . . . . 534.4.1 BCDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.4.2 Joint Gauss-Newton method (JGN) . . . . . . . . . . . . . . . . . . 555 Applications - synthetic examples . . . . . . . . . . . . . . . . . . . . . . . 575.1 Seawater intrusion - state estimation with ADMM . . . . . . . . . . . . . . 575.1.1 Parametrization for the synthetic scenarios . . . . . . . . . . . . . . 595.1.2 Joint inversion with ADMM . . . . . . . . . . . . . . . . . . . . . . . 625.1.3 Joint inversion by ADMM with inexact GW parameters . . . . . . . 665.1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 715.2 Solute tracer test - Parameter and state estimation without an explicit petro-physical relationship . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.2.1 Parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.2.2 Example 1 and 2: Estimation for initial state c0 . . . . . . . . . . . 775.2.3 Example 3: Estimation for hydraulic conductivity K . . . . . . . . . 805.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.3 Joint inversion with electromagnetic loop survey . . . . . . . . . . . . . . . 905.3.1 Joint inversion with ADMM . . . . . . . . . . . . . . . . . . . . . . . 915.3.2 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 936 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95viBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99viiList of TablesTable 5.1 Parametrization for the seawater intrusion Cases 1 and 2. . . . . . . . . . . . 60Table 5.2 Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 for the twodifferent cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Table 5.3 Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 for Gauss-Newton method with the substitution of σ; 1. The weights were fixed throughoutthe minimization based on the initial data misfits; 2. δf was adjusted as δf =φD,eφD,fat each iteration; 3. δe was adjusted as δe =φD,fφD,eat each iteration. GNinversion run up to 10 iterations. . . . . . . . . . . . . . . . . . . . . . . . . 66Table 5.4 Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 whendifferent from true GW parameters are used in the ADMM inversion or coupledapproach. Test 1 - change in permeability field. Test 2 - 70% reduction inpermeability and dispersion values. . . . . . . . . . . . . . . . . . . . . . . . 70Table 5.5 Computational cost overview for the forward model (f.m.) runs and actual runtime of the inversion; all results displayed were executed on a desktop withIntel(R) Core(TM)i7-2860 QM processor and 16GB RAM. Parameters for linesearch procedure were the same for all methods, amount of iteration counts isdisplayed in the brackets next to each method. . . . . . . . . . . . . . . . . . 72Table 5.6 Parametrization for the solute tracer Examples 1, 2 and 3. . . . . . . . . 75Table 5.7 Minimization parameters for the inverse problem . . . . . . . . . . . . . 76Table 5.8 Example 1: Summary of the errors for the solute distribution and electricalconductivity estimates, the data misfits and the cross-gradient coupling terms(between the estimate and the true field) by either separate inversion or jointinversion with cross-gradient product. . . . . . . . . . . . . . . . . . . . . . . 78viiiTable 5.9 Example 1 and 2: Summary of the relative errors (joint inversion vs GW (orDC) data only) for the solute distribution and electrical conductivity estimatesbased on the results from 10 different GW survey designs. . . . . . . . . . . . 78Table 5.10 Example 2: Summary of the errors for the solute distribution and electricalconductivity estimates, the data misfits and the cross-gradient coupling terms(between the estimate and the true field) by either separate inversion or jointinversion with cross-gradient product. . . . . . . . . . . . . . . . . . . . . . . 80Table 5.11 Example 3: Estimation of hydraulic conductivity with KL expansionparametrization: Summary of the errors for the solute distribution, electricalconductivity and hydraulic conductivity estimates (and also the estimate for µ= mean(K), the true value was mu = −1.8) by either separate inversion or jointinversion with cross-gradient product. A final cross-gradient product between cnand σ is displayed too. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Table 5.12 Example 3: Summary of the relative errors (joint inversion vs GW data only)for the solute distribution at cn, hydraulic head h and logarithm of errors forconductivity K based the results from 10 different GW survey designs. . . . . . 84Table 5.13 Errors for the initial and final solute content at time t0 and t1 for the two differentCases 1 and 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91ixList of FiguresFigure 1.1 A simple dipole - dipole survey with one receiver pair and one source pair, from[73]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.2 Coupled hydrogeophysical problem: each problem can be treated separately(blue and orange loop) or the petrophysical relationship can be added and bothmodels updated in some coupled fashion (yellow loop). . . . . . . . . . . . . 6Figure 2.1 A single cell in the staggered grid of the discretized groundwater model, hy-draulic head h and solute concentration c live in the cell centers, fluxes arecomputed at the cell faces. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Figure 2.2 Visualization of the Semi-Lagrangian scheme in 2D: Interpolation of solute massxt+∆t which was advected from the cell center at xt to surrounding cells. Wetrack the particles (cell centers) forward along the flowlines. . . . . . . . . . . 19Figure 2.3 Test against analytical solution for advective transport: Semi-Lagrangianmethod (left) vs Lax-Friedrich implicit scheme (right). Simulations were runfor three different time steps, from top to bottom, t = 15, 40 and 80 days, cor-responding to Courant number 0.33, 1 and 1.33. The solid line is the modeledbreak through curve and the dashed one is the analytic solution. . . . . . . . 25Figure 2.4 The boundary conditions setup for a Henry’s benchmark problem, seawaterboundary is on right, freshwater inflow on the left. . . . . . . . . . . . . . . 26Figure 2.5 Contour plots for two cases of testing the Henry’s benchmark problem represent-ing the steady state solutions for solute mass fraction. Left: the dimensionlessparameters a = 0.33, b = 0.01; right: a = 0.33, b = 0.1, which correspond to thesolutions in [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.6 DC survey scheme for a single source and receiver. . . . . . . . . . . . . . . 28xFigure 2.7 Staggered discretization for the TDEM model, properties such as permeability µand electrical conductivity σ are at the cell centers, magnetic field b is computedat the cell face and electric field e on edges. . . . . . . . . . . . . . . . . . . 30Figure 3.1 The reference model for SWI inverse problem. The resulting conductivity σrefhas high solute mass fraction along the seaward boundary on the left, as ex-pected in real conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 4.1 Comparison of the JGN and BCDM scheme. Left: The scheme for a single iter-ation of the Joint Gauss - Newton minimization, the two data misfit terms ΦD,regularization terms ΦR and a similarity term ΦS are computed for the initialestimates of m and σ, both models are then updated via a Gauss-Newton iter-ation (Eq.(4.34)). Right: BCDM minimization starts with the Gauss-Newtondescent for the groundwater model m only (Eq.(4.32)) involving also the sim-ilarity term, the new GW model mk+1 serves as an input for the geophysicalGauss-Newton descent (Eq.(4.31)), leading to an update for the geophysicalmodel, σk+1, the next BCDM iteration starts with mk+1 and σk+1. . . . . . . 56Figure 5.1 Effect of increased pumping on the propagation of the SWI front. Source:A.D. Wenner, P.E. Jacobsen, L.K. Morgan - Understanding seawater intrusions[poster]. 2013 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Figure 5.2 The geological layers for the heterogeneous case based on Kidd2 site in FraserRiver delta [71], schema for kz field. . . . . . . . . . . . . . . . . . . . . . . 60Figure 5.3 Seawater intrusion scenario: Upper left and right: True models for initial andfinal solute distribution for Case 1; Bottom left and right: True models forinitial and final solute distribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and0.75 are displayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 5.4 Data collection. Left: Experimental design for the DC survey: Dark blue pointsrepresent the electrodes on the surface placed on a regular grid, saltwater iscoming in from the west boundary. Right: GW wells locations, 8 wells in totalare placed along two transects at x = 16 m and x = 24 m distance, later labeledas b and c, other experimental designs for GW data are plotted in Figure 5.10. 62xiFigure 5.5 The error and residual decrease during the ADMM descents: Green trianglesrepresent the scaled error for the final solute fraction ωf , orange stars corre-spond to updated rk values, where rk = σf (k) − p(ωf (k)). The GW wells inthis case were placed along x = 16 m and x = 24 m . . . . . . . . . . . . . . 63Figure 5.6 ADMM estimates for Case 1: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full con-tour lines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5(green) and ω = 0.75 (red). . . . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 5.7 ADMM estimates for Case 2: Contour profiles at x = 10 and x = 20. Thedashed lines are the estimates from the joint inversion (ADMM), and the fullcontour lines are for true locations corresponding to ω = 0.25 (blue), ω = 0.5(green) and ω = 0.75(red). . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 5.8 ADMM estimates. Upper left and right: Estimates for initial and final solutedistribution for Case 1. Bottom left and right: Estimates for initial and final so-lute distribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and 0.75 are displayed.The true models are plotted in Figure 5.3 . . . . . . . . . . . . . . . . . . . 65Figure 5.9 The relative error decrease for five different GW sampling designs; the dottedline is for the ω0 relative error, the full line for the ωf relative error decrease.In all cases the error decrease slows down with further iterations. The plottedresults are based on different transects of wells plotted in Figure 5.10. . . . . . 67Figure 5.10 The transects along x and y axis in plan view. For the homogeneous Case1 thesampling depths were z = [3, 7, 11], and for the heterogeneous Case 2 z = [5, 9].each design had two vertical or horizontal transects. . . . . . . . . . . . . . . 67Figure 5.11 Altering the GW parameters: Errors decrease for ωf estimate and residual rkwhen correct and altered GW parameters (D, k) are used in ADMM. The errorsfor ωf for the coupled approach are plotted as well with a dashed line, it is justa single value. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 5.12 Altering the GW parameters: The true permeability field in x direction kx inthe second test. For solving the inverse problem a fixed value kx = 4.4−11 wasused instead, while the true field is shown in this figure. The true ky and kzwere altered similarly when simulating GW data. . . . . . . . . . . . . . . . 69Figure 5.13 Altering the GW parameters: Errors decrease for ωf estimate and residual rkwhen correct and altered GW permeability is used in ADMM. The errors forωf for the coupled approach are plotted with a dashed line. . . . . . . . . . . 69xiiFigure 5.14 Solute tracer test. Left: A scheme for the GW survey design with twotransect of wells, applied in Example 1 and 2. Right: A scheme for DCsurvey with receivers on surface only, in the background is the electricalconductivity model applied in Example 2 which spatially correlated moreconductive field than background. . . . . . . . . . . . . . . . . . . . . . 75Figure 5.15 Different GW surveys for solute tracer test. Ex 1,2 - des1: Two transects ofwells at x = 14 and 18m applied in Example 1 and 2. Ex 1,2 - grid: 4 x 5 gridof wells from which 8 locations was chosen randomly for each run. Ex3-des1: 3x 3 grid of wells locations for results presented in Table 5.11. Ex3-grid: Regulargrid of well from which 10 wells were randomly chosen, corresponding resultssummarized in Table 5.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 5.16 Solute tracer test. The true models for the initial solute concentrationin Example 1 (left) and Example 2 (right), isosurfaces displayed at c0 =0.2, 0.4, 0.6 and 0.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Figure 5.17 Example 1 and Example 2: Bar plots for the error estimates by different meth-ods, the actual values are in Table 5.8, Table 5.10 respectively. . . . . . . . . 78Figure 5.18 Example 1: The true models vs the estimates by GW data inversion andstructure-coupled inversion, contours at z = 8 and y = 18 for final and ini-tial solute cn, c0 are plotted. . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 5.19 Example 2: The true models vs the estimates by GW data inversion andstructure-coupled inversion, contours at z = 8 and y = 18 for final and ini-tial solute cn, c0 are plotted. . . . . . . . . . . . . . . . . . . . . . . . . . . 81Figure 5.20 Example 2: The estimates for 10 different experimental designs, altering loca-tion of 8 wells from Ex 1,2 - grid design. . . . . . . . . . . . . . . . . . . . . 82Figure 5.21 Example 3 : Bar plot for the errors by different methods, the actual values arein Table 5.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 5.22 Example 3: The true models vs GW data and structure-coupled inversion es-timates for the final solute concentration, contours at z = 8 and y = 16 areplotted and hydraulic head at y = 24. . . . . . . . . . . . . . . . . . . . . . 85Figure 5.23 Example 3: The estimates for 10 experimental designs by different methods,altering the location of 10 wells in Ex3-grid design. . . . . . . . . . . . . . . 86Figure 5.24 Example 3: The true field and estimates of hydraulic conductivity K by GWdata only and JGN and BCDM. . . . . . . . . . . . . . . . . . . . . . . . . 87xiiiFigure 5.25 The experimental design for EM survey: Loop (transmitter) is placed in thecenter on surface, the receiver coils are placed on surface in an uniform grid. . 90Figure 5.26 The error and residual decrease during the ADMM descents: Green trianglesrepresent the scaled error for the final solute fraction ωf , orange stars corre-spond to updated rk values, where rk = σf (k) − p(ωf (k)). The GW wells inthis case were placed along x = 16 m and x = 24 m . . . . . . . . . . . . . . 91Figure 5.27 ADMM estimates for Case 1: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full con-tour lines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5(green) and ω = 0.75 (red). . . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 5.28 ADMM estimates for Case 2: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full con-tour lines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5(green) and ω = 0.75 (red). . . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 5.29 ADMM estimates. Upper left and right: Estimates for initial and final solutedistribution for Case 1. Bottom left and right: Estimates for initial and final so-lute distribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and 0.75 are displayed.The true models are plotted in Figure 5.3 . . . . . . . . . . . . . . . . . . . 93Figure 5.30 The relative error decrease for different GW sampling designs; the dotted lineis for the ω0 relative error, the full line for the ωf relative error decrease. Theplotted results are based on different transects of wells plotted in Figure 5.10. . 94xivAcknowledgementsI would like to start off by thanking my advisor, Eldad Haber, for giving me the opportunityto join the PhD program at UBC, helping me with funding, allowing me to research thetopics I enjoyed most, and guiding me in a research field which was new to both of us. Thefact that I could pursue my PhD studies in Vancouver had a huge impact not only on myprofessional, but also, personal development.I also want to acknowledge Roger Beckie who always found time for me and gave meadvice when I asked. I enjoyed being a teaching assistant in his Groundwater hydrologyand Contaminant transport classes because it gave me an opportunity to explore a role inthe academic world other than research.I want to thank all GIF members for all our valuable discussions and, just as importantly,for creating a pleasant environment to work in. And a special thanks to Mike Mitchell forreviewing this thesis.Last but not the least, I would like to thank to my family for their encouragementduring this long term commitment, but also to all of my “Vancouver families” in capoeira,dance and sailing for the great times outside of work.Special thanks to my mum for her limitless psychical support, to Matei for encouragingme to apply to UBC and always having a good advice when I was frustrated, and to myfriends - Sonja, Romain, Sarka, Daja, Radka, Olda, Ludek, Jonas, Gata, Luz, Dean, Nickand many others for their cheering and support.xvChapter 1IntroductionMany different processes in groundwater hydrology can be labeled as solute transport,comprising different scales and both natural and human affected processes. At a small scalewe can think of fertilizers moving through the root zone, a contaminant plume observedin a few wells can already affect an area of watershed size, or seawater intrusions thatcan affect large regions along the coastal zones. While all these examples evolve fromdifferent applications, they have a few important commonalities. In all these cases weneed a groundwater flow model to determine the flow field in the subsurface and a solutetransport model, to describe, and ideally also to predict the fate of the solute.Theoretically, we do have a groundwater flow and solute transport model specificallytailored for each different process; however, conceptualizing from the real environmentto a mathematical model brings a lots of uncertainty. The soil and porous media arerarely homogeneous, and processes such as dispersion and diffusion are hard to upscalefrom the real scale heterogeneity to the actual model resolution. Furthermore, we need tomake decisions regarding the boundary and initial conditions, which are sometimes fairlystraightforward, but in more complex situations the particular choice can substantiallyimpact the model outcomes.From the numerical point of view, the modeling challenges arise in solving the governingequations, which are a system of two partial differential equations for flow and solute mass1balance:∂ (φρ)∂t+∇ · (ρq) = ρQgw (1.1)∂(φ ρω)∂t+∇ · (φ ρD∇ω)−∇ · (ρω q) = Qsω, (1.2)where ρ is the fluid density, ω is the solute mass fraction, φ porosity, D represents the hydro-dynamic dispersion, t time and Qsω and Qgw the external fluxes of solute and groundwater.These two equations are coupled via groundwater flux q;q = −kµ(∇p+ ρg∇z), (1.3)where k is permeability of the porous media, p fluid pressure, µ fluid viscozity, z is downwarddepth coordinate and g gravity acceleration. Hence, if the fluid density does not change, wecan basically decouple the system and the groundwater flow is essentially governed by thepressure gradient, ∇p, only. For flow of variable density though, the governing equationsare strongly coupled and weakly nonlinear. Any 3D simulation is therefore computationallycostly especially when considering the inverse problem.To overcome the challenge of conceptualization, i.e. simplifying the real world intoa representative physical model, we need to have enough data to properly calibrate thegroundwater model, and hence reduce the uncertainty regarding parameters and initialor boundary conditions. However, data availability is usually somewhat limited for a fewreasons.Obtaining measurements via groundwater sampling can be costly and difficult. Directgroundwater data are usually obtained from boreholes or wells, which even with multiplescreening intervals provide only a few data points information for groundwater states, suchas solute concentration or hydraulic head. Moreover, for variable density flow models addi-tional density data has to be obtained next to standard head measurements. The densityinformation can come from concentration data (respectively electrical conductivity of watersamples), but this is prone to error due to the difference between resident and flow concen-tration (i.e. the collected samples might not represent well the surrounding environment[15]). Besides, parameters like dispersivities or permeability can not be measured directlyat the field scale and we thus have to rely on indirect data and laboratory samples.In many cases, where the complexity of the solute transport process is high, the avail-2able groundwater data is too scarce for a successful calibration and good predictabilityof the underlying models. In these cases, some of the properties such as fluid conductiv-ity, saturation or porosity, can be well captured by geophysical methods sensitive to theseproperties. Indeed, as we show in this work, geophysical methods provide an attractiveoption to add more information for groundwater model calibration, and improve the solutetransport monitoring and predictions.1.1 Geophysical data in hydrologyRecent advances in geophysical instrumentation as well as computational software havemade many geophysical techniques more accessible for environmentally oriented surveys[53]. In groundwater hydrology geophysical methods have already became standard tech-niques to map the bedrock topography, the geology, estimate hydraulic conductivities bysalt tracer tests, and for monitoring seawater intrusions. A new discipline, hydrogeophysics,developed focusing on the use of geophysical measurements for mapping subsurface features,estimating properties and monitoring processes that are important to hydrological studies[84].For example, Direct Current (DC) resistivity methods found a broad application acrossdifferent areas of hydrology. In DC experiments a current is injected into the ground cre-ating an electrical potential, which is then measured using pairs of electrodes - receiversplaced on the surface or in boreholes (Figure 1.1). The DC method was successfully appliedin many hydrogeophysical studies for solute transport, such as tracer tests and seawaterintrusions ([47], [32],[45],[44],[20],[51] or [49]) but also to investigate unsaturated zone trans-port, such as soil salinization [52], or hyporheic flow transport [98]. The conductivity ofsaltwater is a few orders higher than that of freshwater in the pore spaces. Hence, saltwaterintrusion increases the soil bulk electrical conductivity, which makes it an excellent targetfor electromagnetic and electrical resistivity methods. DC surveys were also successfullyused to identify clay lenses, which have higher conductivity compared to other geologicalmaterials, and to determine variable moisture content [12, 73].Electromagnetic methods (EM) are based on the electromagnetic induction effect, andare also sensitive to changes in subsurface conductivity; however, the different physicalprocess of charging conductive bodies in the subsurface can make these methods betterfor recovering the soil bulk conductivity. Depending on the type of survey one can collect1D data, so called EM soundings, with possible large depth investigation or target shallow3Figure 1.1: A simple dipole - dipole survey with one receiver pair and one source pair, from[73].depths with shorter offsets and obtain two and 3D data. In airborne surveys, a helicoptercarrying transmitter and receiver coils can easily map large areas. EM applications ingroundwater hydrology are similar to those of DC methods, but for large scale studies EMsurveys are usually preferred due to easier, though more expensive, data acquisition. Forexample in [77] they collected airborne electromagnetic data in the Red River basin ofTexas to estimate the volume of water intruded by saline water, capturing an area of 34km2. Other examples of seawater intrusion monitoring via EM methods can be found in[72],[93] or [14].Ground penetrating radar is another non-invasive geophysical technique used by hy-drologists, which is sensitive to changes in dielectric permittivity. The reflection depth cangive information about location of lower density plumes [22] or geological structures suchas fluvial deposits [24], which can then help in estimating the hydraulic conductivity ofthe aquifer. Seismic or nuclear magnetic resonance methods have also been used to studyhydrological problems [21, 40, 96].The main advantage of geophysical methods lies in their non-invasive character andability to map a large subsurface area. As opposed to groundwater sampling methods thatgive only point wise observations, geophysical methods can provide extensive 3D data sets4in a relatively short time for a large and heterogeneous area. Many geophysical methodsare available, which together with different experimental designs offer large to small scaleinvestigation of the subsurface. However, the main caveat of geophysical methods is that thecollected data, regardless the chosen method, are indirect. Therefore, an inverse problemhas to be solved to recover the subsurface geophysical property field. This inverse problem isusually ill-posed leading to non-unique solutions dependent on some “a priori” information.Furthermore, the groundwater property of interest is linked to the geophysical property viasome petrophysical relationship, which is based on empirical laws and has to be calibratedfor specific field conditions. The parametrization for this relationship then brings anothersource of uncertainty when “incorporating” the geophysical information into groundwaterestimates. Geophysical data thus come with a price of solving an inverse problem anddetermining the link between the geophysical and groundwater variables.1.2 Hydrogeophysics: different approaches to couple themodalitiesCombining the geophysical and hydrological data increases the amount of data about theflow in subsurface and hopefully can help to overcome some of the caveats of each specific,individual type of data. Ideally, we want to keep the advantages of both types of dataand improve the estimates of groundwater states and parameters. We should keep inmind though that both the geophysical and groundwater inverse problems are ill-posed;with limited data we want to recover the 3D subsurface properties. Hence, both inverseproblems need some form of regularization (a-priori information), that would steer theestimates towards physically realistic ones.As Figure 1.2 shows, each inverse problem can be treated separately, we can add someregularization and solve it for one type of data only. The geophysical estimate can thenbe converted to a groundwater variable via some petrophysical relationship allowing geo-physical and hydrogeological computations to stay independent. The advantage of such anuncoupled framework is that the hydrogeological and geophysical models run independently.However, this also means that the a-priori information from hydrology is not integrated intothe geophysical inversion. Since geophysical inverse problems are ill-posed and require aregularization term, or a prior (if Bayesian methods are used), ignoring hydrogeology datacan severely deteriorate the quality of the geophysical inversion estimates.Another option is to include the coupling term which represents the link between the5Figure 1.2: Coupled hydrogeophysical problem: each problem can be treated separately(blue and orange loop) or the petrophysical relationship can be added and both modelsupdated in some coupled fashion (yellow loop).groundwater and geophysical variables in the inversion. This is usually referred to as acoupled approach. The coupling term can be the previously mentioned petrophysical rela-tionship, which is usually empirical (e.g. Archie’s law which relates soil bulk conductivityand fluid conductivity) or it can be set more loosely by only imposing a structure similarity,between the groundwater and geophysical models.For the actual coupling, if we have an explicit relationship, the geophysical variable(estimate) can be transformed into a groundwater variable or the other way around. Thus,beside evaluating the similarity between the two models, one can also use one of the es-timates as a reference model, and possibly also as an initial guess for solving the otherinverse problem. Another option is “data extraction” [8], where we use the petrophysicalrelationship to convert the geophysical estimates into groundwater states and create someadditional data for groundwater model inversion. For the actual inversion, some approachesinvert only one type of data, however, both data misfits are evaluated with respect to theupdated estimate ([45],[5]).The different ways that the two types of data are combined and solved leads to vari-ous approaches in hydrogeophysical studies and therefore terminology vary among authors.However, in all coupled approaches, the geophysical and groundwater models are linkedtogether during the inversion. The hydrological state estimates are then guaranteed to bephysically realistic and fewer a-priori assumptions need to be entered into the geophys-ical inversion implicitly. Due to this fact, coupled approaches were repeatedly found togive better results ([44], [49] or [47]). However, coupled approaches are computationally6more intensive and possible errors introduced by the coupling term, may map from thegroundwater to geophysical estimate and vice versa and need to be considered.The joint approach can be considered as a specific type of a coupled approach, whereboth data misfits are minimized together. Since both models behave differently in termsof convergence, weighting the two different data misfits next to regularization and couplingterm in one objective function can be a challenging task. Each of the problems can bedifficult to solve by itself due to governing equations, which are partial differential equations,and solving them in a coupled manner can create even more computational complexity.Moreover, the forward and inverse simulations for the geophysical and groundwater modelsare commonly solved by separate codes, possibly “black-box”, which further complicatesthe joint inversion. Therefore, in most literature we find applications of coupled approachesrather than an actual joint minimization.1.3 Coupled hydrogeophysical inversion: The state of artBefore reviewing some of the work on coupled hydrogeophysical inversions, we should recallthat the groundwater (GW) forward problem alone may require solving a system of twostrongly coupled and nonlinear partial differential equations, and thus any 3D simulation iscostly when solving the inverse problem, where multiple forward simulations are necessary.To decrease the amount of estimated parameters and make the inverse problem solvable,usually some a-priori information based on geology is considered to reduce the size of theproblem. One option is to apply some geostatistical constraints when estimating the GWparameters [46, 51, 80]. The actual minimization is then often directed by some generalinverse software such as PEST [25] or UCODE [78] where the sensitivities are derived by aperturbation approach.In the last decade, numerous studies have investigated the potential of coupled inversionand the resulting improvement in parameter estimation when compared to the inversion ofa single data type. For example in [47] the authors tested the coupled versus uncoupledapproaches for an infiltration event on synthetic data by joining HYDRUS1D model andERT data simulations to estimate the soil properties. The authors confirmed the benefits ofthe coupled approach, although the results show that model estimates are highly dependenton how well the hydrologic model represents the real system.The study of [44] also examined different approaches for coupling groundwater and geo-physical data. In their synthetic study a simple cross-sectional groundwater model was7coupled with TDEM (1D soundings) and a real field dataset coupled a Modflow ground-water model with 1D ERT data. Both sequential and joint approaches were tested. Inthe sequential approach the geophysical inversion was independent, but the subsequentgroundwater model inversion used the geophysical estimates, while in the joint approachthe geophysical and groundwater models were updated simultaneously. Geometric andpetrophysical constraints were used in this study to couple the GW model and geophysicaldata by relating a particular layer resistivity to its hydraulic conductivity. For the jointapproach, the authors employed a petrophysical constraint as an additional regularizationterm. They observed that the joint approach increases the computational time comparedto the sequential approach, with no significant improvement of the groundwater parameterestimates in the case of a real dataset.The first joint hydrogeophysical inversion to estimate the hydraulic conductivity appearsin [79] for a synthetic 2D experiment; the authors combined the temporal moments ofelectric potential perturbations with hydraulic head data. Later, the same approach wastested using a real data set from a sandbox laboratory experiment [80]. The fully coupledapproach managed to successfully recover the solute transport and preferential pathways,the hydraulic conductivity fields were smoothed out but contained the main heterogeneousfeatures.In the study of Irving and Singha [49] the authors also jointly inverted geophysicalERT and groundwater concentration data for a saline tracer test; however they introduceda stochastic approach for this type of application. The pilot points technique was usedto generate possible hydraulic conductivity fields in 2D and a Markov chain Monte Carlosampling algorithm was implemented within the joint inversion. In [51] they expanded onthis work by adding self potential data. The synthetic case results for their 2D profileshowed improved estimates when combined data were used compared to an inversion ofeach type of data independently. The authors noted the necessity to parallelize the code toreduce computation time, if applied to 3D problems.For seawater intrusions (SWI), with variable density flow modeling, the situation isusually more complex in terms of time scale and heterogeneity compared to solute tracerexperiments. However, similar ideas were implemented in the calibration of underlyingsolute transport models, leading to a high diversity of approaches when solving the hydro-geophysical inverse problem. In most cases groundwater and geophysical data are invertedseparately and the coupling happens via a petrophysical relationship.In the study by [8] the authors investigated a sequential approach to estimate hydraulic8conductivity and dispersivity parameters, which they tested on two synthetic benchmarkproblems representing pumping experiments in a coastal aquifer in 2D. The ERT derivedconductivities were transformed via Archie’s law to salt mass fraction. These estimateswere then filtered using a cumulative sensitivity based on the squared Jacobian and servedas data in the hydrogeological parameter estimation next to the “collected” salt massfraction data. The inversion was performed with PEST using a gradient based method.The geophysical ERT data did improve the SWI model calibration.The calibration of a seawater intrusion model using geophysical data can be also foundin [45]. In this work the GW model is used to interpret the data and guide the geophysicalinversion. Saltwater concentrations based on the 2D SEAWAT model were converted withArchie’s law to electrical resistivity and the forward modeled response of a TDEM soundingwere calculated and compared with observed data. The TDEM data sets (1D soundings)were collected in Santa Cruz County, California. In total, six parameters were estimatedfor both the groundwater and TDEM model using the PEST optimization system. Theauthors concluded that the coupled approach provided a significant improvement in spatialresolution which would be hard to obtain with standard geophysical inversions.Large scale hydrogeophysical problems are the focus of the recent work by [18]. Theauthors improved the inversion framework of iTOUGH2 to enable parallel computing andmerge it with a parallel geophysical simulator for electromagnetic data, creating a generalframework they can apply to a wide range of processes in multiphase flow and solutetransport. The sensitivities with respect to parameters were evaluated by a perturbationapproach, and the high computational burden of this approach was balanced by the factthat the perturbed model simulations could be run independently.Compared to other synthetic studies and applications, we set up both models in thesame computational environment in order to test different minimization schemes using ajoint approach, which is one of the goals of this thesis. For the joint approach we proposea couple of methods that can split the minimization into two subproblems and thus easethe computational burden of the coupled problem.1.3.1 Structure - coupled inversionMost of coupled inversions and even inversion of geophysical data alone, involve some kind ofpetrophysical relationships that relate the two different physical properties. The petrophys-ical relationship represents an explicit way to impose a similar structure on two different9models. In environments with well mapped geology and well calibrated parameters it isfeasible to generate a petrophysical map that links groundwater and geophysical parame-ters or states, yet such information is generally not available. Some authors question thevalidity of the petrophysical relationships (i.e. knowledge of the porosity distribution). Insome cases researchers are able calibrate it in the field, but often this petrophysical relationis accepted just by applying some fixed empirical parameters for either synthetic or fielddatasets. If incorrectly assigned this direct link may cause significant errors in the solutecontent estimates or enforce the coincidence of the two models, when there is none. Thisnaturally leads to the idea of replacing the fixed petrophysical relationship by a more loosestructural constraint.In the structure-coupled inversion the two or more different models are not linked by anexplicit petrophysical relationship but instead, some similarity is assumed, such that thetwo different models exhibit changes in the same regions. This similarity is imposed by astructural constraint in the minimization procedure leading to model estimates of differentproperties that share the same “structure”.The idea of a structure-coupled inversion to image the subsurface using geophysicalmethods has been studied over two decades [34, 39, 64]. Typically, indirect data are collectedby multiple methods, but due to the different physical phenomena of each method, themodel estimates often vary [35]. Moreover, each set of data might have a very differentresolution and be contaminated by different measurement errors. Still, inverting thesedatasets in some coupled fashion is desired if all data are related to the same geologicalfeature. Examples of structure-coupled inversions with geophysical data can be found in[34, 48, 63, 69]. Joint total variation is particularly popular method for this type of jointgeophysical inversion.Similarly, when coupling the hydrological and geophysical inversions, different ap-proaches exist that alleviate the need of the fixed petrophysical relationship. For examplein the work of [13] the authors assumed that the parameter of interest was mainly affectedby geological facies and used a level set method to minimize the coupled objective function.The same method could be expanded for joint inversion with groundwater and geophysicaldata using the same level set parametrization (geological facies distribution). In [54], theauthors replaced the petrophysical relationship by an assumption of a strong correlationbetween the fluid and soil bulk electrical conductivity. In particular, for their syntheticstudy with a solute tracer test, they used ERT time lapse data together with groundwaterhead and fluid conductivity data to determine hydraulic conductivity.10To our knowledge, the first work applying a structural constraint in hydrogeophysi-cal inverse problem is the study of Lochbuhler et al.[65]. The authors combine ground-penetrating radar data with hydraulic tomography or tracer mean arrival times to estimatehydraulic conductivity. As a structural constraint they applied a cross-gradient field prod-uct. In this thesis we further expand on the ideas of [65], but implement a full 3D inversionand jointly invert geophysical and hydrological data that are not directly dependent on theproperty of interest.1.4 Thesis objectives and research contributionsThe motivation for this thesis was to improve estimates for modeling of solute transportprocesses by combining hydrogeological and geophysical data. This objective was conqueredin few different steps, which can be summarized as follows:• Developping the groundwater and geophysical models in conjunction withanalytically derived data sensitivitiesJoint approaches were reported to improve the groundwater estimates. However,the forward and inverse simulations for the geophysical and groundwater models arecommonly solved separately to take advantage of each inversion code and forwardmodel. In the joint approach both data misfits are minimized together, which of coursebrings computational complexity by implementing the coupled framework, solving alarger problem and weighting the two different data misfits next to regularization andcoupling term in the objective function.In order to overcome some of these issues, we developed both models in the same com-putational environment, so that we can test different minimization schemes for a jointinversion of hydrological and geophysical data. The common framework also simplifies thecoupling implementation between the physically different models. We derived the sensitiv-ities of data with respect to parameters or states analytically, using a closed form formula.Unlike the perturbation approach, the analytical derivation of sensitivities significantly re-duces the computation burden, and thus enables us to solve problems in 3D. We testedboth models against benchmark problems or by tests using artificial sources.• Propose computationally efficient inversion schemeFor the sake of simplicity and testing the approaches themselves, the hydrogeophysi-cal inversions are often either of a small size or applying 2D data only. However, in11the field we can rather expect large domain and dataset, especially in case of geo-physical data. In this thesis we therefore proposed a computationally efficient scheme,which can handle large datasets (e.g. discretization with more than 200 cells in eachdirection)Specifically, we introduced a version of the alternating direction method of multipliers(ADMM), which allows us to efficiently split the objective function of the coupled inverseproblem. ADMM was introduced in early 70’s and has recently gained popularity for manyinverse problems. It is a natural choice for multiphysics problems [11], due to its strongconvergence properties [36]. ADMM is suitable for cases where the constraint (here thepetrophysical relationship) can be considered as exact, in practice having low uncertaintiesusually suffices. A successful hydrological application can be found in Wohlberg et al. [99],where ADMM was applied to solve the inverse problem to estimate the piece-wise smoothhydraulic conductivity fields. Many more applications can be found in machine learning orstatistical modeling.In addition to ADMM we also examined other minimization techniques such as theBlock Coordinate descent method (BCDM), which also splits the minimization onto twosubproblems and a standard Gauss-Newton method, that updates both models at eachiteration. Both Block coordinate descent and ADMM have a major advantage over theGauss-Newton method. First, by splitting up the objective function into geophysical andgroundwater part, we do not need to weight two different data misfits with respect toregularization and coupling term at each iteration. There are almost no criteria for settingtwo parameters during the iterative minimization process when each data misfit has adifferent speed of convergence. Second, the multiplication of Jacobians and data misfit canhave a very different computational cost for each model and hence make each iteration veryunbalanced and difficult to parallelize.The ADMM method was tested on a seawater intrusion scenario, we jointly inverted thesolute mass fraction data together with either DC resistivity or TDEM data to estimatethe groundwater states.• Structure-coupled inversionAnother limitation of the joint hydrogeophysical inversion is the knowledge of thepetrophysical relationship. Despite its abundant use, this empirical relationship ishard to properly parametrize in field conditions and adds uncertainty to the coupled12inversion estimates. In this thesis we further explored the potential of joint hydrogeo-physical inversion when the parameters for petrophysical relationships are unknownand implemented structural constraints such as the cross-gradient field product andjoint total variation instead, so that the hydrological and geophysical data can be stillinverted in a joint approach.The coupling via structural constraints usually evaluates the gradient fields of eachproperty rather than the actual values and thus encourage the two models to change in thesame location. Therefore, it is a good alternative to a petrophysical link, as it does notenforce an exact linear mapping between soil bulk conductivity and solute content whenthere is none.We tested the structure coupled inversion on a synthetic test of solute tracer to estimatethe initial conditions for solute content and hydraulic conductivity field. Groundwatersolute concentration and hydraulic head data were jointly inverted with DC resistivitydata.In the following chapter we elaborate on both groundwater and geophysical modelsused in this work, their governing equations, our chosen discretization and our approachfor the solution of the forward problem. In chapter 3 we show a general scheme for theminimization of a single type of data, which has crucial importance when building up moreadvanced joint inversion schemes. The actual coupled inversion is detailed in chapter 4,where we also expand on coupling options such as petrophysical relationships and structuresimilarity constraints, and introduce the ADMM method together with other methods forjoint inversion. Applications of all proposed methods are documented in chapter 5; section5.1 is devoted to the seawater intrusion scenario (results of this experiment were publishedin [91]) and section 5.2 shows a structure-coupled inversion for a solute tracer test (themanuscript is under review). In the last chapter we summarize our findings.13Chapter 2Forward modelingIn this chapter we present the models that we developed to test different minimizationschemes in the hydrogeophysical inversion. For solute transport we have two models thatcan simulate solute tracer tests and seawater intrusions scenarios. Both of them are basedon solute and fluid mass balance equations. For the geophysical modeling, direct currentresistivity and time domain electromagnetic data were simulated. For each model we discussthe governing equations, applied assumptions, discretization and solution of the equationsand present the results of some model verification tests.2.1 Groundwater modelModels of solute transport can range from simplistic to fairly complex representations of thereal world. The choice often depends on the application, required accuracy and availabledata for model parametrization. In solute tracer tests or plume monitoring, the time scaleis often short, and density variation effects do not need to be considered. However, inapplications such as seawater intrusions a full variable density flow model is desirable toimprove the accuracy of forecasts.In our study, we developed two groundwater models, both consist of the fluid and solutemass balance equation, where each one is represented by a partial differential equation in14time and space:∂ (φρ)∂t+∇ · (ρq) = ρQgw (2.1)∂(φ ρω)∂t+∇ · (ρω q)−∇ · (φ ρD∇ω) = Qsω (2.2)The system above is a pressure - solute mass fraction formulation [57], where ρ is thefluid density, ω is the solute mass fraction, φ is porosity, D represents the hydrodynamicdispersion, t is time and Qsω and Qgw are the external fluxes of solute and groundwaterrespectively. The two governing equations are coupled via groundwater flux q (sometimesreferred to as specific discharge), which, taking the general formulation for Darcy’s law [7],can be written asq = −kµ(∇p+ ρg∇z), (2.3)where p is the hydraulic pressure, k is the permeability of the porous media, µ is the fluidviscosity, g is the gravitational acceleration and z is an upward coordinate direction.Depending on the problem we want to model, we need to add further assumptions orsimplify the system above. If the fluid density changes across the model domain due tothe solute content dynamic, the system becomes nonlinear and one needs to complete itwith state equations which define the dependence of density ρ and other porous mediaparameters on the solute mass fraction. However, if the solute dynamic does not effect thedensity we can simplify the system by taking the density ρ out of the brackets. Furthermore,for steady state conditions, the velocity term q does not change and the two PDEs can besolved separately.In the next two sections we describe the specifics of a simple solute tracer model and avariable density flow model, and elaborate on the methods applied to solve the discretizedsystems of equations.2.1.1 Modeling of solute tracer movementFor the solute tracer model we assume fluid of a uniform density, therefore, we can rearrangeequation (2.3) (Darcy’s law), and use the formulation with hydraulic head and conductivityinstead. Hydraulic head is essentially the sum of pressure and elevation head (potential)15and in unconfined aquifers it corresponds to the water level in the well:h =pρg+ z. (2.4)After taking the density term ρg out of the brackets in Eq.(2.3) we obtain:q = −kρgµ(∇ pρg+∇z) = −K ∇h, (2.5)where K is hydraulic conductivity, which is equal to kρgµ and has length per time units,[L/t], unlike the permeability k which has [L2] units and is related to the architecture ofthe porous media only.In this application we assumed steady state conditions. This assumption can not begenerally applied, however, if we look closer at the time dependent term and consider thefluid to be compressible, we can apply chain rule and rewrite the time dependent term as:∂(φρ)∂t= φ∂ρ(p)∂t+ ρ∂φ(p)∂t= φ∂ρ∂p∂p∂t+ ρ∂φ∂p∂p∂t=(φ∂ρ∂p+ ρ∂φ∂p)∂p∂t. (2.6)Letting α be the coefficient of (rock) matrix compressibility, α = ∂φ∂p , and β, the coefficientof fluid compressibility, β = 1ρ∂ρ∂p ; the specific storativity is then Ss = ρg(φβ + α) [6] andwe can write the time dependent term for hydraulic head as(φ∂ρ∂p+ ρ∂φ∂p)∂p∂t= ρS∂h∂t=1gSs∂p∂t= ρSs∂h∂t. (2.7)Ss, specific storativity, is the volume released (or added) from or to a unit volume of storagedue to a unit decline (increase) in hydraulic head. The units are [L−1].Since the solute tracer model is for flow of uniform density, after dropping out thedensity term and substituting (2.5) in (2.1) we obtain:Ss∂h∂t−∇ · (K ∇h) = Qgw (2.8)In confined aquifers specific storativity has fairly low values of order 10−3 or lower [92]. Wecan therefore neglect the time dependent term Ss∂h∂t , when∣∣Ss ∂h∂t ∣∣ << |Qgw|, i.e. the exter-nal sources representing pumping rates or boundary effects are much larger than changesdue to hydraulic head decline (rise). This will be mostly the case for a short time scale and16confined aquifers with pumping tests.The steady state fluid mass balance equation (2.8) thus becomes−∇ ·K∇h = Qgw (2.9)Appropriate boundary conditions need to be defined, here we use Neumann and Dirichletboundary conditionsh = hBC at ΓD,−→n ·K∇h = 0 at ΓN . (2.10)hBC is the fixed head along the boundary ΓD, and a no-flux boundary is set for ΓN . Thegroundwater linear velocity v, which is the velocity of solute moving in the subsurface, isdefined asv =1φq = − 1φK∇h. (2.11)Now we can rearrange the solute mass balance equation (2.2), and take the densityterm out of brackets, which leads to a standard advection-diffusion equation for soluteconcentration c (replacing the ρω term). Assuming, that porosity φ does not change intime or vanish we obtain the following time dependent partial differential equation for cwith groundwater velocity v:∂c∂t+∇ · (c v)−∇ · (D∇c) = Qsw. (2.12)This is a standard advection-diffusion equation. Even though that in aquifer systemsthe groundwater velocity is fairly small, for the scenarios we consider (and model) in thiswork, it is the dominant process compared to dispersion.Discretization and solution of the solute tracer modelBoth governing equations (2.9) and (2.12) are discretized on staggered grids, with c andh placed in the cell centers and fluxes at the cell faces (see in Figure 2.1). The velocityfield (specific discharge) q is computed at the cell faces to avoid long differences whencomputing the spatial gradient of h and improve accuracy of the discretization scheme.The cell faces values of porous media properties such as dispersion, hydraulic conductivityK or porosity, that are necessary to calculate the velocity v, are obtained in our model byharmonic averaging from the values in cell centers.Discretization of the spacial derivatives is done using a finite volume approach. For17Figure 2.1: A single cell in the staggered grid of the discretized groundwater model, hy-draulic head h and solute concentration c live in the cell centers, fluxes are computedat the cell faces.time integration we use operator splitting, where we first solve the advective part with aSemi-Lagrangian integrator and then parabolic part with an implicit Euler method (fordetails see [37]). The discretized system for the solute tracer groundwater model is then:A(K) h = qgw (2.13)cn+1 − S(v)cn∆t= A(D) cn+1 + qsc, (2.14)where A(K) = −Div Av(K)Grad, A(D) = Div Av(D) GradThe matrices Div and Grad represent the divergence and gradient operators and Av(·)harmonically averages the values of a model property in the x, y and z direction from cellcenters onto cell faces. S(v) represents the Lagrangian push forward operator moving thesolute along the flow lines given by velocity v (discussed in the next paragraph), ∆t is thetime step size and external flux qgw accounts for the pumping rates and also the fixed headboundary effects. Similarly qsc represent the sink/source term for the solute.The solution of fluid mass balance is straightforward; the hydraulic head h from Eq.2.13 is simply obtained ash = A(K)−1qgw, (2.15)which can be solved directly or by Conjugate Gradient method for large scale problems.The solute mass balance equation (2.14) is solved by operator splitting. The choice ofa Semi-Lagrangian method here for the advection part is motivated by the ability to takean arbitrary large time step which is especially important for inverse problems, where the18velocity field may not be known. Standard explicit techniques may require changing thetime step during the inversion process which could lead to instabilities in the computationof the gradient.Our Semi-Lagrangian method belongs to the family of modified methods of charac-teristic (MMOC) first introduced in [26]. The main advantages are in the alleviation ofthe Courant number restriction due to the Lagrangian advection step [16] and mass con-servation. The Eulerian - Lagrangian scheme is also effective in overcoming numericaldispersion for advection dominated problems [89]. Similar approaches are taken for ex-ample by codes such as the MOC3D model for solute transport [58], or later MT3DMSfor variable density flow in connection with Modflow [60, 62]. In the context of reviewon Eulerian-Lagrangian localized adjoint methods (ELLAM)[85] our approach is a finitedifference Eulerian-Lagrangian type, where we do not solve the solute transport equationusing an integral equation but instead use a projection matrix.In our implementation, particles are placed at the cell centers at each time step andtracked forward. The solute content of each transported particle is interpolated to itsneighboring cells yielding the solution at the next time step (see in Figure 2.2). TheLagrangian push forward operator S(v) is thus formed at each time step by coefficientsthat interpolate the transported particles (cell centers) onto the neighboring cells.Figure 2.2: Visualization of the Semi-Lagrangian scheme in 2D: Interpolation of solute massxt+∆t which was advected from the cell center at xt to surrounding cells. We trackthe particles (cell centers) forward along the flowlines.19Afterward, we can add the dispersion step using Euler’s method and equation (2.14)can be solved recursively in time as:cn+1 = (I −∆tA(D))−1 (S(v)cn + ∆tqsc) . (2.16)Insight for Lagrangian methodsThe basic idea behind Semi-Lagrangian methods can be well explained using a 1D advection equation witha fixed velocity:∂c∂t+ v∂c∂x= 0The analytical solution is c(x, t) = c0(x− vt). Suppose we have a discretization over a uniform mesh wherexi = iδx and tn = nδt, the Lagrangian methods use this analytic solution to find c at time tn+1 as:c(xi, tn+1) = c0(xi − v (n+ 1)δt) = c0(xi − vδt− vnδt) = c(xi − vδt, tn)Similarly forc(xi, tn) = c0(xi − v nδt) = c0(xi + vδt− v nδt) = c(xi + vδt, tn+1)For the Lagrangian methods we can therefore obtain an exact solution; however, at each time step we wouldhave to change the mesh. The Semi - Lagrangian methods approximate instead the value based on theneighboring points using some interpolation technique. For 1D linear interpolation we can simply writecn+1i as:cn+1i = αcni−k−1 + (1− α)cni−kwherevδtδx= k + α, k =⌊vδtδx⌋.For a 3D case we can proceed similarly, using trilinear interpolation, where the coefficients α and 1−α willcreate the matrix S and we can compute c at the next time step ascn+1 = S cn.For a 3D domain with n cells in each direction, S is n3 x n3 matrix containing the interpolation coefficientscalculated using some interpolation procedure (e.g. linear interpolation leads to first order accuracy). Theoverall accuracy of the advection step depends on the chosen interpolation method when deriving thecoefficients for the S matrix.For the advection case a straightforward connection exist to the so called upwind scheme which isdefined as:cn+1i = cni − v δtδx{(cni+1 − cni ), v < 0(cni − cni−1), v ≥ 0with the stability region |v| δtδx≤ 1. Thus for |v| δtδx≤ 1 the upwind scheme is the same as Semi-Lagrangianscheme above.202.1.2 Variable density flow (VDF) modelThe modeling of variable density flow (VDF) has various applications including contaminanttransport, saline lakes hydrology, brine migration or nuclear waste disposal monitoring andany processes, where the fluid density changes with solute concentration or temperature.The physics of VDF can follow two very different processes. First, if the density of the fluiddecreases with depth, or significantly denser fluid enters the system above lighter fluid, itmay lead to unstable mixing processes (so called fingering) and for such cases the governingequations for VDF become highly nonlinear. On the other hand if the density increaseswith depth the process of mixing and flow is stable. This is the case for most groundwatersystems where the density differs laterally but usually increases with depth [95]. In somecase, for example salt tracer tests, the density changes are small enough not to effect theflow field and a standard solute transport model can be used instead.VDF models are used in this study for seawater intrusion simulations, where the differ-ence between freshwater and saltwater density is relatively small; however, the changes inthe flow field due to solute dynamic still require the governing equations to be solved in acoupled manner.For the variable density flow the situation time dependent term is a bit different thanin the case of solute tracer model. We can again assume that porosity changes only withpressure, but the density now changes with both solute mass fraction and pressure, followingthe chain rule the term ∂(φρ)∂t can be written as∂(φρ)∂t= φ∂ρ(p, ω)∂t+ ρ∂φ(p)∂t= φ∂ρ∂ω∂ω∂t+ φ∂ρ∂p∂p∂t+ ρ∂φ∂p∂p∂t(2.17)= φ∂ρ∂ω∂ω∂t+(φ∂ρ∂p+ ρ∂φ∂p)∂p∂t. (2.18)Similarly as in section 2.1.1 we can rewrite the fluid mass balance equation (2.1) usingspecific mass storativity Smp = ρ(α+ φβ) =1gS as [6]:φ∂ρ∂ω∂ω∂t+ Smp∂p∂t+∇ · (ρq) = ρQgw (2.19)For confined aquifers we can expect∣∣∣Smp ∂p∂t ∣∣∣ << |ρQgw|. The first term is equal to φργ ∂ω∂t ,where γ is introduced in the state equation bellow (Eq.(2.22)) and for saltwater-freshwaterdensity difference it is usually bellow 0.025. Therefore in many cases∣∣φργ ∂ω∂t ∣∣ will be much21smaller than |ρQgw| and the steady state assumption can be justified.The governing equations for variable density flow, with steady state assumption forpressure, are then∇ · (ρq) = ρQgw, with q(ρ) = −kµ(∇p+ ρg∇z), (2.20)∂(φ ρω)∂t+∇ · (ρω q)−∇ · (φ ρD∇ω) = Qsω. (2.21)We complete the system with the linearized state equation between fluid density ρ andsolute mass fraction ω:ρ = ρF (1 + γω) with γ =ρS − ρFρF, (2.22)where ρF is the freshwater and ρS the saltwater density.Additional state equations can be considered for modeling VDF, the effect of pressureand temperature changes on density or the effect of solute content on viscosity µ. However,for our application and under isothermal conditions, the viscosity changes can be neglectedand the solute content dynamic will have a much higher effect on density than pressurechanges. Another simplification of our model involves keeping the hydrodynamic disper-sion tensor D fixed, i.e. not dependent on Darcy’s velocity q, which slightly reduces thecomputational burden, when solving the coupled governing PDEs iteratively. Still, using afull variable density flow model, as for example in [57], would not change the approachesdescribed in the following chapters regarding solving the coupled inverse problem.The two PDEs (2.20) and (2.21) are coupled via the velocity term q which is dependenton density ρ, respectively the solute mass fraction ω. The forward VDF model is thereforenonlinear and both equations need to be solved together. We can again apply a finitevolume approach for spatial derivatives; averaging, gradient and divergence operators todiscretize the fluid mass balance equation, which becomes:Div RKM (Grad pn + (Avρn) g) = ρn Qgw (2.23)where RKM = diag(k Av(µρ)−1), g = −g Grad z.We use the Hadamard product a b for the element wise product of two vectors and Avmatrix operator for harmonic averaging from cell centers onto cell faces. The ad-script n22stands for the corresponding hydrological state at the n − th time step. For now, assumethat we know the density ρ at the nth time step, the pressure pn can be then solved directlypn = (−Div RKMGrad)−1 (Div RKM (Avρn) g − ρn Qgw) , (2.24)by using Cholesky decomposition for matrix −Div RKMGrad.Given the pressure pn we can compute the groundwater linear velocity vn at the cellfaces:vn = −diag(Av(µφk))−1(Grad pn + (Avρn) g) . (2.25)For the solute mass balance equation (2.21) we can again apply operator splitting;however, it is not the solute concentration c this time but the mass of the particle ρω whichis tracked forward during the advective step. The discretized form is then(ρ ω)n+1 − Sn(ρ ω)n∆t= Div A(D,ρn+1) Grad ωn+1 + Qsω (2.26)with A(D,ρn+1) = diag(Av1D ρ(ωn+1))−1.Here again Sn represent the Semi-Lagrangian operator, containing the interpolation weightsbased on velocity field vn. First, we integrate the advection step explicitly, obtaining sometemporary solute mass (ρ ω)∗:(ρ ω)∗ = Sn(ρ ω)n (2.27)Using the state equation (2.22), we can solve a local quadratic equation for ω∗, noting thatonly one root of the equation makes physical sense. Having the solution ω∗, we can nowintegrate the diffusion part, starting from ω∗ with an implicit Euler method as:(ρ ω)n+1 = (ρ ω)∗ + ∆tDiv A(D,ρn+1) Grad ωn+1 + Qsω (2.28)The implicit diffusion step is nonlinear and is solved using a Picard iteration [2], updatingthe velocity field vn and corresponding ω∗. For variable density flow with seawater, thedensity difference between the two fluids is fairly small, decreasing the nonlinearity of thecoupled system and the number of Picard iterations (see the Algorithm 1). Even thoughthe time step can be large given the stability of the Semi - Lagrangian method, care must23be taken with respect to the coupling with the flow equation [81]. Too large of a step wouldlead to a weak coupling and possibly erroneous calculations.Algorithm 1 Picard schemeInitialize ωn+11 = ωn,ρn+11 = ρnFor k = 1 to itermax• Calculate of pressure field: pn+1k = f(ωn+1k );.• Update velocity field: vn+1k+1 = − kµ(∇pn+1k + ρn+1k g∇z))• Advection step: (ρ ω)∗ = Sn(ρ ω)n• Dispersion step: (ρ ω)n+1k+1 = (ρ ω)∗ + ∆tDiv A(D,ρ∗) Grad ωn+1k+1 + Qsω.• Update pressure field given the new solute mass fraction: pn+1k+1 = f(ωn+1k+1)• Check the error: ∥∥pn+1k+1 − pn+1k ∥∥ < tol break else Step12.1.3 Testing the GW modelsTo verify the developed groundwater models, we first tested numerically the discretizationerrors. We computed the first and second derivatives for a known solution, which was thencompared with the model result for a decreasing grid size. From the hydrogeological pointof view we tested our model simulations against a classical benchmark problem for VDFand known solutions for a 1D advection diffusion equation in case of solute transport.If we fix the velocity field for the system of equations in the solute tracer model, ana-lytical solutions exist for the solute mass balance equation (2.12), either for advection onlyor with some specific initial conditions also for advection-diffusion process. Therefore, wecan compare our model results with this analytical solution and also evaluate the accu-racy of Semi-Lagrangian method against the Lax-Friedrich scheme. In a simple experimentwe set the hydraulic head boundary conditions to be Dirichlet for the inflow and outflow,implying a uniform velocity field along the x axis. With such a setup the solute with afixed concentration equal to 1 at the inflow moves with a steady velocity from left to right,and without any dispersion we would have a sharp interface. We solved this essential 1Dadvection problem with our 3D model for different time steps using both Semi-Lagrangianand Lax-Friedrich methods. The results (profiles of the concentration profiles) are shownin Figure 2.3, confirming the stability and lower numerical diffusion of the Semi-Lagrangianmethod.To test the VDF model, we run our code for the Henry problem [43] which is a clas-sical benchmark problem for VDF in 2D representing a simplified seawater intrusion case.24Figure 2.3: Test against analytical solution for advective transport: Semi-Lagrangianmethod (left) vs Lax-Friedrich implicit scheme (right). Simulations were run for threedifferent time steps, from top to bottom, t = 15, 40 and 80 days, corresponding toCourant number 0.33, 1 and 1.33. The solid line is the modeled break through curveand the dashed one is the analytic solution.25The boundary conditions are set as impermeable on the top and bottom of the domain,the inland boundary is set as fixed freshwater inflow, and the seaside boundary is fixedhydrostatic pressure based on the seawater interface (see in Figure 2.4). We used the di-mensionless parameters a = 0.3214 and b = 0.1 as in [1], and obtained similar results as intheir study for a diffusive case of the Henry problem, i.e. when only molecular diffusion isconsider with a fixed value for D and transverse and longitudinal dispersivities being zero(see in Figure 2.5).Figure 2.4: The boundary conditions setup for a Henry’s benchmark problem, seawaterboundary is on right, freshwater inflow on the left.Figure 2.5: Contour plots for two cases of testing the Henry’s benchmark problem repre-senting the steady state solutions for solute mass fraction. Left: the dimensionlessparameters a = 0.33, b = 0.01; right: a = 0.33, b = 0.1, which correspond to thesolutions in [1]262.2 Geophysical modeling2.2.1 Direct current resistivity surveyIn order to estimate the physical properties of the porous media, we choose direct current(DC) resistivity, also referred to as Electrical Resistance Tomography (ERT). DC resistivityis sensitive to the electrical conductivity of the media. Since the conductivity of saltwater isa few orders of magnitude larger than that of freshwater, the DC resistivity method is oftena natural choice, where the fluid conductivity changes as a result of the solute dynamic.Moreover, DC has several advantages when compared to other electromagnetic methods:relatively easy and inexpensive data acquisition and minimal computational effort.In DC experiments, a current is injected into the ground creating an electrical potentialdistribution, which is then measured using pairs of electrodes placed on or under the surface(Figure 2.6). Multiple receiver and source electrodes can be utilized in a surface survey,creating many options for the experimental design scheme, with the goal of capturing thevariability of the 3D subsurface conductivity distribution.To model this process, we used the steady state form of Maxwell’s equations:∇ · (−σ∇ϕ) = I(δ(r− rs+)− δ(r− rs−)) (2.29)n · ∇ϕ = 0 on Γncwhere σ represents the media’s electrical conductivity, ϕ is the electric potential, I is thecurrent source, δ is a Dirac delta function and rs+,s− stands for the location of positiveand negative current electrodes. The boundary conditions were set as no flux across theboundaries, Γnc. When solving the forward problem, the electrical conductivity model, σ,is known and potentials everywhere can be calculated using a finite volume approach on a3D grid. Since the discretized DC equation is essentially the same type of equation (Poissonequation) as the flow equation (2.13), we follow the same procedure to solve it as it hasalready been described in section 2.1.1. The cell-centered finite volume approach for thediscretization of the problem leads to a linear system of equations:A(σ)u = q (2.30)where A(σ) = −Div diag(1(Avσ−1))Grad.27The electrical conductivity σ is averaged harmonically from the cell centers onto cell faces(via matrix operator Av), q is the source term and u represents the potentials. The for-ward model solves for the potential everywhere given a conductivity model. Using a dataprojection matrix Qe, we obtain observed datade = Qeu (2.31)measured at the receivers.Figure 2.6: DC survey scheme for a single source and receiver.2.2.2 Time domain electromagnetic survey (TDEM)All electromagnetic surveys employ the induction effect that creates magnetic field due tochanges in electric field and vice versa. The way how the changes are initiated and thedata are recorded give different type of surveys, here though, we focus just on one specifictype that leads to a particular set of governing equations and following discretization ofthe electromagnetic process.In our experiments, a direct current runs through the transmitter loop to create asteady-state magnetic field. The current is then abruptly shutoff, causing a rapid changein magnetic field, which induces so called eddy currents that flow in electrically conductivebodies in the subsurface. These eddy currents will generate secondary magnetic fields thatcan be measured by receivers at the earth surface. Similarly as DC resistivity, TDEMsurvey is thus sensitive to changes in electrical conductivity of the subsurface.Maxwell’s equations for the electric field e and magnetic density b, where we already28substitute in the constitutive relations, can be written as:∇× e+ ∂b∂t= 0 (2.32)∇× µ−1b− σe− ∂e∂t= sr(t) (2.33)and BC:~n× b = 0, b(0, x) = b0, e(0, x) = 0 (2.34)where e is the electric field, b is the magnetic flux, µ represents the magnetic permeability,σ is the electrical conductivity, electrical permittivity and sr is the source current. Thefirst equation is often referred to as Maxwell’s Faraday equation, and the second one asAmpere’s law.We can apply the quasi static approximation here due to the initial shutoff of the source.Since the permittivity is very low and we can assume small changes of electric field e, theterm ∂e∂t can be dropped out and for t > 0 (after the current shutoff) we can simplify theequation (2.33) to:∇× µ−1b− σe = 0 (2.35)Discretization and solving the forward modelAfter elimination of b in Eq.(2.35) by substituting e from Eq.(2.32), we can write:∇× µ−1∇× e− σ∂e∂t= 0, (2.36)where the only unknown is the electric field e.Using a staggered grid, the disretized form of the governing equations (2.32) and (2.36)in time is then:C en + α(bn − bn−1) = 0 (2.37)C> M(µ)−1C en + αM(σ)en = αM(σ)en−1 (2.38)where M(µ) = diag(Afµ−1) M(σ) = diag (Avσ) .The electric field is placed on the edges and the magnetic density occurs on the cell faces(see in Figure 2.7). C is the curl matrix operator, α(i) is 1δt and Ae,Av correspond toaveraging matrix operators from cell centers onto edges and cell faces respectively.29Figure 2.7: Staggered discretization for the TDEM model, properties such as permeabilityµ and electrical conductivity σ are at the cell centers, magnetic field b is computed atthe cell face and electric field e on edges.In order to proceed with the forward calculation in time, we need to know the initialfield e0 and magnetic flux b0. The initial magnetic field is calculated using the Biot-Savartlaw and integrating over the current ellipse (loop). Substituting into Eq.(2.35) we obtainthe initial electric field ase0 = M(σ)−1(C>M(µ)b0).Once we have the initial fields, we need to solve the time dependent Eq.(2.38) for e, at eachtime step we thus solve the following system for ei+1:(C> M(µ) C + αiM(σ))ei+1 = αiM(σ)eiThis system is solved with the Conjugate Gradient method, for which we implemented theA-phi preconditioner [37]. The magnetic density at i + 1the time step is obtained directlyfrom Eq.(2.37) asbi+1 = bi − 1αiC>ei+1.The “measured data”, changes in magnetic field are then:di = Q C ei+1, (2.39)where the projection matrix Q interpolates the data from the cell faces onto receivers30locations.Both geophysical models were tested using artificial sources. These were analyticallyderived based on known conductivity fields, which could be then compared against modelsolutions.31Chapter 3Inverse problem for a single modelThe focus of this work is on the joint inversion of two different data sets. However, before westart to explore the beauty of joint inversion, it is worthwhile to at least briefly summarizesolving the inverse problem for a single data set since it is a building block for any coupledminimization scheme. We restrict the two following sections to a general solution of aninverse problem for a single model, but leave some space to discuss specifics related toeach, groundwater and geophysical, type of data. In the last section of this chapter wederive the data sensitivities, which are necessary if we want to proceed with any type ofNewton method for a single or coupled inversion.3.1 Ill-posedness and regularizationUsually the data we collect provide some indirect, less often direct, information about themodel we are trying to recover. One could try to estimate this model by considering onlythe data alone and minimize the norm of the data residuals:minmΦ =12‖Q F(m)− dobs‖2 , (3.1)where F(m) is some forward model operator, e.g. mapping the model m to hydrologicalstates or geophysical potentials (and most likely involving some inverse matrix computa-tion), matrix Q projects the simulated groundwater or geophysical states onto the datalocations and the vector dobs represents the observed data for some groundwater or geo-physical states which are dependent on the model m that we are trying to recover. Letsassume a linear case, where we can write Am = QF m, the least square solution of Eq. (3.1)32would then give mˆ [66]:mˆ = (A>A)−1A>dobs (3.2)If the data are coming from some physical experiment, are scarce and contain measurementerrors, the minimization will rarely lead to success. Moreover, the forward operator F(m)can behave well for a simple groundwater model but it is rarely well behaved for geophysicalmodeling. With indirect data for m, the matrix (A>A) may not be invertible. Even if weobtain direct GW data, for example solute concentration, the Q operator has nonzeroentries only for few data location points in the observation wells, making the recovery of afull 3D model difficult (i.e. nonunique solution).Both groundwater and geophysical data provide only a limited “picture” of the sub-surface, whereas we want to recover a full 3D model, whether it is solute concentration orsoil bulk conductivity. Therefore, we need to add some a-priori information to limit thesolutions to physically realistic estimates, which are hopefully close to the true ones.Thiscan be achieved using regularization. The inverse problems are usually regularized by as-suming that the model (estimate) should be close to some reference model, which can besimply an initial guess or derived from geological information. Clearly, the way we chooseour a-priori information and incorporate it in the regularization will have a huge impact onthe final estimate.If we set R(·) to be some regularization operator and Σ to be the covariance matrixbased on measurement errors and the norm ‖r‖2Σ−1 = r>Σ−1r, we can obtain the maximum-a-posteriori estimate of the model m by minimizing the following functional:minmΦ =12‖Q F(m)− dobs‖2Σ−1 + αR(m), (3.3)The regularization parameter, α, should be chosen in such a way that the inversion is guidedby the data misfit term. In other words, we prefer the data to be the steering mechanismin the minimization rather than the a-priori information, which may be biased or incorrect.The parameter α can be set by a trial and error or by Cross-Validation [37].Since the model to be estimated in our application is based on differing fluid conductivitywe do not expect to have sharp boundaries due to dispersion in porous media. Therefore,for our application, L2 quadratic regularization turned out to be a good choice, as it favors33smooth solutions. After discretization the quadratic R(m) givesR(m) =12‖L(m−mref )‖2 , (3.4)where L is the gradient operator, and mref is the reference model. The matrix operatorL consists of gradient fields in x, y and z directions; L = [αxGradx αyGrady αzGradz]>,and each of the components can be proportionally decreased or increased if one expectsmore smoothness in a particular direction. For example, in the case of seawater intrusionswith solute coming from the seaward boundary, the variation along this boundary canbe penalized more than other directions. The selection of the reference model is up topersonal preference, in geophysics a uniform model is a popular option, unless we havemore information about the geology. The reference model can be also updated during theiterative minimization process.An important thing to note here is that often we do not estimate the conductivityor solute concentration directly but rather estimate some function of m. For examplef(m) = exp(m) ensures that the resulting field attain only positive values. In the case ofseawater intrusion modeling the tanh(x) function turned out to be useful to parametrizethe solute concentration. The electrical conductivity, resp. solute fraction, was modeled asσ = 1− tanh(m),where mref varies only along x coordinate (flow direction) and is bounded so that σ staysattains only positive values. In Figure 3.1 we show a slice through the reference modelalong the x axis, the reference model mref increases gradually along the x direction, whilethe resulting σref mimics the seawater intrusion expected profile (i.e. solute coming fromthe left seaward boundary).3.2 Gauss-Newton methodLet us recall the objective function in equation (3.3), after substituting the quadratic reg-ularization for m we have:Φ(m) =12‖Q F(m)− dobs‖2Σ−1 +α2‖L(m−mref )‖2 ,34Figure 3.1: The reference model for SWI inverse problem. The resulting conductivity σrefhas high solute mass fraction along the seaward boundary on the left, as expected inreal conditions.which can be minimized by a Gauss-Newton method. The gradient of Φ(m) is thengm = J>W>W (F(m)− dobs) + α L>L (m−mref) , (3.5)where J represents the sensitivity of data with respect to model m, J = ∂F (m)∂m , and W =diag(σˆ−1)Q. The Hessian can be approximated byHm = J>W>WJ + α L>L. (3.6)The computation of the sensitivities J is the subject of the following section.35At each iteration we update the estimate mk by a step in the direction:sk = −H−1m gm, (3.7)which is solved using a preconditioned conjugate gradient (CG) method. The advantage ofa CG method is that it does not require computing the gradient or Hessian explicitly, onlymatrix vector products need to be computed. The CG method thus eases the computationalburden especially for large scale problems or when solving the problem on a fine scale. Forfurther acceleration of the iterative process it is wise to choose a good preconditioner; fora single type of data we used an approximation of the inverse of ∂2R(m)∂m2, the regularizationpart in Hessian. The step size µ is determined afterward using an Armijo line search [75]giving the update for m as:mk+1 = mk + µsk. (3.8)The number of iterations differs for groundwater and geophysical data, as the problemshave different convergence behavior.3.3 SensitivitiesTo apply any type of Newton method, we need to know the sensitivities, or so calledJacobians, of collected data with respect to parameters or states of interest. In this sectionwe show how to derive them analytically based on the discretized systems of equations,using a closed form formula. Such analytically derived sensitivities enable us to speed upthe inversion process and ease the computational burden.For the solute transport model (the system of governing equations was presented in thesection 2.1.1) we need to compute the groundwater data sensitivities: ∂c∂m and∂h∂m , wherem represents either hydraulic conductivity K or initial solute concentration c0.The derivation of sensitivities of hydraulic heads h w.r.t. hydraulic conductivity K isquite straightforward, the fluid mass balance in Eq.(2.13) can be rewritten as F(K, h) = 0and differentiating w.r.t. K we obtain∂F∂K+∂F∂h∂h∂K= 0.36Rearranging, we have∂h∂K= −(∂F∂h)−1 ∂F∂K= −A(K)−1Div diag (Grad h)Av. (3.9)The sensitivities of solute concentration data w.r.t. initial solute c0 can be solvedrecursively from the discretized equation (2.14):∂cn+1∂c0= (I −∆tA(D))−1 S∂cn∂c0. (3.10)Jacobian w.r.t to the initial conditions, c0, is therefore solved as a by-product when per-forming the forward calculations. The adjoint can be solved similarly by a backward timestepping process.It is important to note that by using Krylov subspace methods (CG in particular) forsolving the direction update in Gauss-Newton method (Eq.(3.7)), we do not need to storethe explicit sensitivity matrix J and can instead work only with matrix-vector products: Jvand adjoint JTw.In particular, for a given vector v and sensitivity ∂cn∂c0, at each kth step we first calculateS∂ck−1∂c0v, then we solve a linear system (I −∆tA(D)) x = S∂cn∂c0v for x at each time step,and the solution x at the last time step gives Jv. Similarly for the adjoint product J>w,we start from the last time step and solve system (I −∆tA(D))> y = w, and then multiplythe solution y by S>. After all backward time stepping we obtain J>w. Therefore, for thesensitivity computation we only need to solve a linear system corresponding to the size ofvector cn.For the variable density flow model (section 2.1.2) the system of governing equations isnonlinear, which complicates the sensitivity derivation. One way to look at it is that duringthe forward simulation we solve the system of two partial differential equations, and at theend of each time step the pressure pn is given by the solute distribution and boundaryconditions. The groundwater velocity vn can be therefore expressed as a function of ωnand pn only and the system thus reduces to the second equation for solute transport (2.26).For the sake of simplicity, lets proceed using the solute mass fraction formulation alone,the time stepping process in equation (2.26) can be then written as(I−∆tA(D))ωn − Sn−1ωn−1 −∆t qn−1BC,ex = 0. (3.11)37Due to the operator splitting approach, we can explicitly derive the sensitivity of the massωn+1 at each time step with respect to ωn, the solute fraction at the previous time step.This is a two step calculation, step 1 being the advection part and step 2 being the diffusionpart. The sensitivity of the final solute fraction can therefore be calculated recursively dur-ing the forward groundwater model run. However, unlike in the solute transport case, theSemi-Lagrangian matrix operator Sn is now dependent on velocity vn, which is dependenton ωn. This extra term makes the derivation more tedious.Let Tn be the sensitivity of Sn(ρω)n with respect to spatial coordinates:Tn =∂(Sn(ρω)fixed,n)∂xi.The sensitivity of the advected mass w.r.t ω0 can be then expressed as∂(ρω)n+1,ad∂ω0= Sn∂(ρω)n∂ω0+ ∆t Tn∂vn∂ω0. (3.12)The sensitivity of ∂vn∂ω0can also be (after some work) expressed as linearly dependent on∂(ρω)n∂ω0since vn is a function of pressure pn and density ρn and thus the term ∂(ρω)n+1,ad∂ω0can be solved recursively.More generally, we can also consider all time steps together:F(ω) =I −∆t A(D)−S1 I −∆t A(D)−S2 I −∆t A(D)...−Sn−1 I −∆t A(D)ω1ω2....ωn−S00....0ω0 −∆tq1BC,exq2BC,ex....qnBC,ex = 0,We can then write the forward model in a compact form asATSω − B0ω0 − q̂ = 0, (3.13)where ATS is the time stepping matrix that is block bidiagonal and B0 is the matrix thatmultiplies ω0. The vector q̂ involves the boundary conditions and sources and vector ω isset here as ω = [ω>1 , . . . ,ω>n ]>.38Differentiating F with respect to ω0 we obtain that∂F(ω,ω0)∂ω∂ω∂ω0+∂F(ω,ω0)∂ω0= 0, (3.14)and therefore∂ω∂ω0= −∂F(ω,ω0)∂ω−1∂F(ω,ω0)∂ω0. (3.15)Using (3.13) we see that sensitivity J is∂ω∂ω0= A−1TSB0. (3.16)Since in our synthetic examples we collect data only for the last time step, we computethe sensitivities recursively during the forward simulation.In a similar way, the sensitivities of magnetic flux changes with respect to the con-ductivity σ are solved recursively during the time stepping process. Differentiating thefirst governing equation (2.37) for the time domain electromagnetic modeling (which is notdirectly dependent on σ) we can write:C∂en∂σ+ α(∂(bn − bn−1)∂σ)= 0. (3.17)Afterward, if we differentiate the second governing equation (2.38) for the electric field andmove some of the terms on the right hand side we obtain:C> M(µ)−1C∂en∂σ+αM(σ)∂en∂σ= α∂M(σ)en−1,fixed∂σ+αM(σ)∂en−1∂σ−α∂M(σ)en,fixed∂σ(3.18)The left side has the same matrix multiplication as the original equation (2.37), except thatnow we are solving for ∂en∂σ and the right hand side is dependent only on σ,∂en−1∂σ and en.Therefore, to compute ∂ei+1∂σ we can just apply the same TDEM forward problem solver,but with a different right hand side. The sensitivity of magnetic flux data with respect toσ is then:∂d∂σ= Q C∂ei+1∂σ. (3.19)The adjoint product is solved recursively with backward time stepping.The sensitivities of the potentials measured by DC resistivity method, ∂ϕ∂σ , follow the39same lines as the sensitivity of hydraulic head with respect to hydraulic conductivity inequation (3.9), since in both cases a Poisson type of equation is discretized. Again, onlymatrix-vector products are calculated.40Chapter 4Coupled inverse problemThe idea of jointly evaluating different data goes back as far as 1975, when Vozoff andJupp [97] jointly inverted two different kinds of geophysical measurements (DC resistivityand magnetotelluric data) by considering a simple three layer model. Inverting jointly hy-drological and geophysical data has much more recent history, first studies were publishedabout a decade ago and focused on unsaturated zone, as for example in [59], where theauthors tested the joint inversion approach on two synthetic examples, combining GPRdata with water saturation data. There are two main reasons for a relatively recent his-tory of joint hydrogeophysical inversion. First, while in geophysics often both types ofdata relate to the same physical property, in hydrogeophysics this is usually not the caseand a relationship between the groundwater and geophysical variable of interest has to beestablished. This gives ample space to create a framework for coupling geophysical andgroundwater data, but the particular coupling can also be potentially the main source oferror when jointly inverting these two datasets. Second, the groundwater and geophysicaldata are usually collected and analysed by two different research communities, with limitedcommunication in between. Here, we bridge this gap by implementing both models withinone computational framework. First, in section 4.1, we describe different approaches tocouple the groundwater and geophysical variables, in section 4.2 we formulate the coupledinverse problem, and finally in the sections 4.3 and 4.4 we show different approaches tosolve the coupled inverse problem.414.1 Coupling the GW and geophysical variablesThe necessary condition to proceed with any coupled inversion is that both the geophysicaland groundwater data, directly or indirectly, relate to the same physical property, whetherit is fluid conductivity, porosity, fluid saturation etc. The actual relation is then commonlydefined by some empirical petrophysical relationships that couples the geophysical andhydrological states or parameters, for example electrical conductivity with solute contentor seismic velocity with porosity, fluid conductivity or geological facies location. Suchpetrophysical relationships have to be either calibrated in the field or can be estimatedfrom the laboratory analysis.An alternative to a petrophysical relationship is applying a structural constraint inthe minimization that has a preference for structure similarity between the two models,however, the two models do not need to coincide. The consistency between the modelsideally results from the data and not from the assumed petrophysical relations.In the next two sections we therefore discuss both approaches to couple the geophysicaland groundwater variables; petrophysical relationships and structure similarity. Both ofthese options are later applied in our synthetic examples.4.1.1 Petrophysical relationshipsPetrophysical relationships represent the simplest and probably the most common way tocouple the two different models. Widely used Archie’s law was introduced already in 1942 [3]and provides a link between electrical conductivity, σb, and fluid conductivity σw. Anotherexample of petrophysical relationship is the link between the hydraulic conductivity andcharge density when measuring the self-potentials [50]. In saturated porous media the fluidhas usually the major effect on the overall soil bulk electrical conductivity. At the sametime, the fluid conductivity is often directly related to the amount of solutes in groundwater,which is our quantity of interest. Therefore, in the context of this work, the Archie’s lawrelating fluid and soil bulk conductivity is the most relevant petrophysical relationship andis the subject of the next paragraph.Coupling the salt mass fraction and electrical conductivityThe soil bulk electrical conductivity of a porous media is dependent on fluid conductivity,rock/soil porosity, permeability, saturation, temperature and also mineral composition [56].In the saturated sediments of coastal aquifers, soil bulk conductivity is usually mainly42affected by the fluid (groundwater) salinity, though if the groundwater composition does notchange then surface conductivity of fluid-grain interfaces, porosity or amount of saturationcan be the main factors influencing the overall soil bulk conductivity. For fully saturatedsediments the Archie’s law can be written as:σb =1ασwφn, (4.1)where σb is the bulk electrical conductivity, σw is the fluid electrical conductivity and φporosity. α and n are empirical parameters related to a geological material, which can becalibrated based on soil core samples, a field survey or recommended values from geologicalcharacterization. In this form of Archie’s law the surface conductivity is neglected, whichis a valid assumption for most of clay free sand and gravel aquifers, however, not forenvironments with a high clay content [56]. In many cases though, the differences inconductivity due to the variation in geological material are negligible compared to theincrease in conductivity due to solute dynamic in groundwater.The fluid conductivity in natural waters increases with the amount of dissolved solidsand ions [42]. Linear relationships can be found for a fixed temperature between totaldissolved solids (TDS) and fluid electrical conductivity. For example the seawater is usuallydominated by sodium chloride ions and we can assume a linear relationship between saltmass fraction and fluid conductivity asσw = b ω + σF , (4.2)where σF is conductivity of freshwater, b is a constant and ω is salt mass fraction. Theelectrical conductivity of seawater is approximately two orders of magnitude higher thanthat of freshwater. After substituting into Archie’s law, we haveσb =1α(bω + σF )φn. (4.3)In this petrophysical relationship one assumes that the bulk conductivity is affected only bythe electrical conductivity of the fluid in the porous matrix, the variations in temperatureare small, and that the surface conductivity of porous material is negligible.Equation (4.3) is used to generate electrical conductivity models based on the salt massfraction from groundwater model simulations in the synthetic examples in section 5.1. Forfield applications it is more common to record the fluid conductivity σw of groundwater43samples, and only later by using a linear relationship to transform it to an actual solutemass fraction or concentration. However, this linear relationship is different from the generalpetrophysical relationship, having a different error and lower uncertainty. Therefore, in thiswork, we prefer to relate the solute fraction (or concentration) and soil bulk conductivity.4.1.2 Structure-based couplingAs it has been pointed out in [70], most petrophysical relationships were originally de-fined for a scale of a representative elementary volume and therefore it is generally hardto determine and properly parametrize them for the field conditions. Another caveat ofapplying the petrophysical relations is the assumption that the properties such as porosityor saturation are homogeneous or known, a fixed parameter setting is then unlikely to bevalid over the entire domain in environments with different geological units.To avoid applying the petrophysical relationship one can assume instead that the twomodels have a common interface, or lithological units [70]; however, for such an assumptionwe need to have a good a-priori knowledge about the geology of the area, which may not bereadily available, especially if we are about to model the processes in 3D. Application of thisapproach can be found in [13]. Another way to alleviate the need of explicit petrophysicallink is by assuming a strong correlation between the groundwater and geophysical variables[54].A more robust structure similarity measure seems like a good alternative to imposing astrict dependence by a petrophysical link, assuming a correlation or defining a-priori somefixed geological model. We consider two models m1 and m2 to have a similar structure ifthey change in the same locations (their transition zones coincide) [38]. The sign of changedoes not need to be the same; however, for the models of soil bulk electrical conductivityand solute content this is the case. The review paper of [35] provides different examples ofapplying the structure similarity measures in multiphysics geophysical imaging, where mostof them involve computation of spatial gradients of the two different models (variables).The gradients can be normalized as in [27], serve to detect edge locations [39] via somemodel curvature function or become an input for a cross-gradient function or joint totalvariation.44The edge detection idea [39] lies in defining a polynomial with some threshold τ1, τ2:S(m) =0 if∣∣∇2m∣∣ < τ1P5(∣∣∇2m∣∣) τ1 ≤ ∣∣∇2m∣∣ ≤ τ21 if τ2 ≤∣∣∇2m∣∣ (4.4)This give smooth operator S(m), fifth degree polynomial, which is also twice differentiableand thus convenient for Gauss-Newton type of minimization.Joint total variation (JTV) of the two gradient fields is becoming popular in multimodalgeophysical inversion and is based on the idea of sparse signal joint recovery. It encouragesthe two gradients to occur at the same location. For multiple geophysical data it has beenapplied in [38] or [28], though the idea originates in processing of the color images [10, 88].For two models m1 and m2 it isΦJTV (m1,m2) =∫ √|∇m1|2 + |∇m2|2dx,where m1 and m2 represent the two, possibly different, models with a similar structure.Computing the cross-gradient field product (CGP) is another convenient option to eval-uate structure similarity between two models:ΦCGP (m1,m2) =∫|∇m1 ×∇m2|2 dx.In geophysics it was introduced by [33] to jointly invert seismic and DC resistivity data andin [65] it was first applied for hydrological and geophysical data inversion.Both JTV and CGP similarity measures can be differentiated with respect to bothmodels and JTV is also convex. One thing to note is that if one of the models is flat, (i.e.having a single value everywhere) then its gradient field is equal to zero. The ΦJTV thenchanges to standard total variation regularization operator for the other model, but in caseof ΦCGP the overall similarity term goes to zero.This is important to keep in mind when evaluating the coupled objective function. Ifthe gradient |∇m1| ≈ 0, the other model m2 will not be penalized for any variations asφCGP goes to zero. Additional L2 smooth regularization can reduce these variations butwill be in conflict with creating a “structure”. Therefore, in [38] it was proposed to add45another term in φCGP to stabilize the problem in areas where |∇m1| ≈ 0, which is∫|∇m2|χθ (|∇m1|) dx, where χθ(t) = 12(tanh(θt) + 1) . (4.5)The parameter θ is chosen so that the value of θ |∇m1| is large enough, where we consider∇m1 to be significant.In this work we tested Cross Gradient product and Joint Total variation as structuresimilarity measures. Both of them are easy to implement and can be also differentiated.4.1.3 Discretization of the coupling termsLets start with the petrophysical relationship, which as a part of the objective function willchange toΦPR = ‖σ(ω)− σ‖2 , (4.6)in other words we will minimize the norm of the difference between the soil bulk con-ductivity, and the soil bulk conductivity given by the solute mass fraction. We need todifferentiate it with respect to both models and thus obtain the gradients with respect toσ and ω. These sensitivities can be easily derived as:∇σ ΦPR = −(σ(ω)− σ), ∇ω ΦPR = ∂σ(ω)∂ω(σ(ω)− σ) = bαφm(σ(ω)− σ). (4.7)And the second derivatives are as follows∇2σ ΦPR = −I, ∇2ω ΦPR = (bαφm)2I, (4.8)with I being the identity matrix.The discretization of the JTV term is quite straightforward as well, we just need to keepin mind that the gradient fields occur at the cell faces and need to be averaged back on thecell centers values. Let Grad be the gradient operator and Af averaging matrix operator(from cell faces onto cell centers), the discretized JTV for the two models m1 and m2 isthen:ΦJTV (m1,m2) = vT√M, (4.9)where M = Af((Grad m1)2 + (Grad m2)2)46and v is a vector corresponding to the cell sizes of the domain and sums up the term underthe square root.We can twice differentiate ΦJTV with respect to m1 and m2, which is necessary to solvethe inverse problem by Gauss-Newton method. The gradient w.r.t. m1 is:∇m1 ΦJTV = GradTdiag (Grad m1) Af T M−12 (4.10)and the Hessian is:∇2m1 ΦJTV = GradTdiag (Af TM−12 ) Grad−GradT diag (Grad m1)Af T M−32 Af diag (Grad m1)Grad (4.11)For practical calculations the first term in (4.11) is semi-positive definite and therefore agood approximation to Hessian.We could also directly discretize the cross-gradient field product (CGP), but that isgenerally numerically unstable due to the long differences that are needed to discretize thecross product. Therefore, we avoid the direct discretization of the cross product, and usethe Lagrange’s identity in a vector form instead and write φCGP as:φCGP (m1,m2) =∫|∇m1 ×∇m2|2 dx =∫|∇m1|2 |∇m2|2 −∣∣(∇m1)T∇m2∣∣2 dx. (4.12)In this case only dot products are used, derivatives along the same directions are multipliedand ΦCGP term can be then discretized asΦCGP (m1,m2) =12vT(Af (Grad m1)2 Af (Grad m2)2 − (Af(Grad m1) Af(Grad m2))2), (4.13)where again v is a vector corresponding to the cell sizes of the domain and Af is averagingmatrix operator (from cell faces onto cell centers).Following the equation (4.13), the gradient w.r.t. m1 for cross-gradient field product is∇m1 ΦCGP = vTdiag (AfGrad m2) Afdiag (Grad m1) Grad−vTdiag (Af(Grad m1) (Grad m2)) Afdiag (Grad m1)Grad, (4.14)47and the Hessian is∇2m1 ΦCGP = GradTdiag (AfGrad m2) Af Tdiag (v) Grad−GradTdiag((Grad m2) AfTdiag (v)diag (Grad m2))Afdiag (Grad m2)Grad. (4.15)Both terms were implemented in our inversion code for the Hessian calculation.4.2 Formulation of the coupled problemUp to now we have described how to discretize and solve both groundwater and geophysicalforward models, how to solve a single inverse problem including necessary regularizationand also two different options for coupling the groundwater and geophysical variables. Weare thus in a good shape to proceed to a coupled problem and solve it with a joint approach.To keep the discussion general, we can assume that we are interested in some GWstate or parameter (model) m which can represent the hydraulic conductivity K or initialconditions for solute content (or also external fluxes, dispersion or other parameters). Thegroundwater data are dependent on model m, with s(·) representing some hydrologicalstate (solute concentration or hydraulic head). The geophysical model to be estimated isthe soil bulk electrical conductivity, σ, that is affected mainly by the fluid conductivity andtherefore is expected to have a similar structure as groundwater solute distribution.Assume now that we have obtained both types of data at a single time, that is we havede = Qeu(σ) + e (4.16)df = Qfs(m) + f (4.17)where Qe and Qf are sparse matrices that project the geophysical states, u(σ), and thehydrological states, s(m), onto their measurement locations, respectively. Let e and fbe vectors with the errors associated with each measurement, which are assumed to beGaussian, independent and identically distributed with covariance matrices Σe and Σf .We now develop a procedure to estimate the groundwater model m and the electricalconductivity σ by the regularized maximum likelihood estimate [30, 37, 67]. This leads to48the following constrained optimization problemminσ,m12‖de − Qeu(σ)‖2Σ−1e + αeR(σ) + δ12‖df − Qfs(m)‖2Σ−1f + αfR(m) (4.18)s.t ΦS(s(m),σ) ≤ τ (4.19)Here R is the regularization operator (discussed in section 3.4) and αe, αf and δ are regu-larization parameters. δ specifically weights the contributions of two different data misfits.The constraint (4.19) here represents either petrophysical relationship (Archie’s law) orsome structure similarity constraint between hydrological state s(m) and soil bulk con-ductivity σ. The choice of τ , representing the error bound for the similarity constraint,depends on the particular similarity coupling, if it is set to zero, we impose an exact simi-larity between groundwater and geophysical model. In the next few sections we will discusspossible ways of minimizing this functional.4.3 Minimization with an explicit petrophysical constraintLets first assume we do know the explicit petrophysical relationship between the two models,σ = η s(m) + σfb = p(m), (4.20)which is based on the Archie’s law introduced in the equation (4.3), and therefore s(m)refers here to solute mass fraction. We lumped a few parameters into η = cαφn, and σfbcorresponds to the conductivity of porous media with freshwater only (or without anysolute) and is equal to σfb =1αφnσF . To recall, the objective function for the coupledinverse problem (equation (4.18)) is then:minσ,m12‖de − Qeu(σ)‖2Σ−1e + αeR(σ) + δ12‖df − Qfs(m)‖2Σ−1f + αfR(m) (4.21)s.t σ − p(m) = 0 (4.22)A number of different approaches may be taken to solve the optimization problem above.One possibility is a direct substitution of p(m) for σ and work with m alone. In such casewe solve the problem with an exact constraint and the objective function (where we omit49the extra regularization term for σ) isΦ(m) = δe12‖de − Qeu (p(m))‖2Σ−1e + δf12‖df − Qfs(m)‖2Σ−1f + R(m). (4.23)This approach has the advantage of solving a smaller problem, however, it complicatesother aspects of the inversion. First, we are required to choose two weights, δe and δf , forthe data misfits with respect to regularization term at each iteration. Despite abundantresearch on choosing one regularization parameter, there are almost no criteria for settingtwo parameters during the iterative minimization process when each data misfit has adifferent speed of convergence. A more detailed discussion can be found in [17], where theauthors used magnetotelluric data together with controlled source EM data. Even thoughboth data sets can be modeled by changes in conductivity, non-trivial weighting was neededin order to jointly invert them. Second, the multiplication of Jacobians and data misfit canhave a very different computational cost for each data misfit, which can make each iterationvery unbalanced and does not favor parallelization. Finally, if the relationships between σand p(m) are inexact, forcing them may lead to inversion artifacts. Therefore, rather thansolving the problem “all at once” we decided to apply the alternating direction method ofmultipliers (ADMM) to minimize this coupled objective function.4.3.1 ADMMThe main advantage of ADMM is that the groundwater and geophysical parts can be solvedseparately, that is, we do not need to weight the two different data misfits in one objectivefunction, but instead we split the minimization into two subproblems. This enables us touse existing inversion methodologies and even software packages only with minor changes.It also takes the advantage of the existing parallelization for an inversion of a single model(discussed in the previous chapter) and therefore substantially increase the efficiency whensolving large scale problems. Furthermore, since the petrophysics is enforced as a constraint,it is not applied exactly throughout the path of optimization which yields additional degreesof freedom to the optimization algorithm.Following the ADMM approach, the augmented Lagrangian for (4.21) isL(σ,m,y) = 12‖de − Qeu(σ)‖2Σ−1e + αeR(σ) +12‖df − Qfs(m)‖2Σ−1f + αfR(m)+y>(σ − p(m)) + ρ2‖σ − p(m)‖2 (4.24)50Here y is a Lagrange multiplier and ρ is a parameter that can be chosen somewhat arbi-trarily. The kth iteration of ADMM is summarized in Algorithm 2.Algorithm 2 ADMM• Approximately minimize the augmented Lagrangian with respect to σ.• Approximately minimize the augmented Lagrangian with respect to m.• Update the Lagrange multiplier yk+1 = yk + ρ(σk+1 − p(mk+1)).We now discuss the solution of each subproblem and show that by using small modifi-cations to existing inversion codes, the ADMM iteration can be carried out efficiently.Geophysical imaging block descentAt each step of minimizing the augmented Lagrangian with respect to σ, we approximatelysolve:minσΦ(σ) =12‖de − Qeu(σ)‖2Σ−1e + αeR(σ) + y>(σ − p(m)) + ρ2‖σ − p(m)‖2 (4.25)The objective function consists of a data misfit and a regularization part, as in usual inverseproblems, and coupling terms involving also m, the groundwater variable. We can proceedusing the Gauss-Newton method to minimize (4.25) with respect to σ. Compared to astandard inverse problem, we need to know also the derivatives of the coupling terms withrespect to σ, however, these are, in this case, straightforward since m is fixed and werealready derived in section 4.1.3.If we set Je =∂u(σ)∂σ and assume for simplicity that Σe = I, the derivative∂Φ∂σ is∂Φ∂σ= JeTQeT (Qeu(σ)− de) + αeLTL(σ − σref ) + y + ρ(σ − p(ω)) = 0 (4.26)and the search direction at each time step is∆σ = − (JeTQeTQe Je + αeLTL + ρI)−1 ∂Φ∂σ. (4.27)The model is then updated by a “soft” Armijo line search [75], σk+1 = σk+µ∆σ, similarlyas for a standard inverse problem described in section 3.2. The parameter µ is adjusted toensure sufficient decrease of Φ, where k refers to kth iteration of the geophysical descent. Wetake a small fixed number of Gauss-Newton steps, and as has been already mentioned in 3.4,51we do not search directly for the electrical conductivity σ model m, where σ = 1−tanh(m).Groundwater model block descentSimilarly, the objective function for the augmented Lagrangian with respect to groundwatermodel m, can be written as:minmΦ(m) =12‖df − Qfs(m)‖2Σ−1f + αfR(m) + y>(σ − p(m)) + ρ2‖σ − p(m)‖2 (4.28)Now the geophysical variable σ is fixed and the derivatives of the coupling terms withrespect to m involve the sensitivity Jf =∂s(m)∂m , which were derived in section 3.3.Again, we assume a quadratic regularization for m , which gives the gradient of Φ(m)as∂Φ∂m= Jf>Qf> (df − Qfs(m)) + αfL>L(m−mref )− ηJf>y + ρηJf>(σ − p(m)) = 0. (4.29)The groundwater model minimization follows the same Gauss-Newton approach as has beendescribed above for geophysical imaging.ADMM stopping criteriaThe ADMM algorithm is stopped when the norm of the residual rk, i.e. the constraintgiven by the petrophysical relationship rk = σ − p(m) , is sufficiently small or the changesin ‖rk‖ between the last few iterations are bellow some threshold value. At this point,we have the capability of fitting both the geophysical and groundwater data such that thegeophysical and groundwater model agree.It has been shown that ADMM has a linear rate of convergence1, where for some appli-cations the desired precision can be reached in a relatively small amount of steps/descents[11]. The penalty term ρ has an effect on the speed of convergence [23, 36]; however, undermild conditions any positive value of ρ will lead to convergence [36], both in terms of theresidual rk → 0 and finding an optimal solution for both models. In our study we set upρ by a trial and error procedure, where ρ ∈ [0.1, 1] provided similar results. In some casesbetter rate of convergence can be achieved by so called Over-relaxed ADMM, or by adding1Note that Gauss-Newton has a linear rate of convergence as well although the constant in Gauss-Newtonmay be better than the ADMM constant.52a scaling parameter for the Lagrangian multiplier, the work of [74] provides actual ratebounds based on the parameter choices.4.4 Minimization without an explicit petrophysicalconstraintAs we already elaborated on in section 4.1.2, the petrophysical relationship used in ADMMminimization is often not known a-priori, its parametrization is uncertain or has to be cali-brated. If we instead apply some structure similarity measure ΦS , we can again change theconstrained formulation of the coupled objective function in Eq.(4.18) to the unconstrainedone and obtain:min Φ(m,σ) =12‖Qeu(σ)− de‖2Σ−1e + αe ‖Le(σ − σref )‖+δ12‖Qfs(m)− df‖2Σ−1f + αf ‖Lg(m−mref )‖+ βΦS(s(m),σ)(4.30)We can apply a Gauss - Newton method and minimize Φ(m,σ) simultaneously for bothmodels, the soil bulk conductivity σ and groundwater model m. However, care has tobe taken with regard to a different convergence speed and data misfit magnitudes of thegroundwater and geophysical models. An alternative is to apply a Block Coordinate DescentMethod (BCDM), which, similarly as ADMM, splits the coupled objective function intotwo subproblems, and treats each separately with a Gauss-Newton method.4.4.1 BCDMBlock Coordinate Descent method minimizes each subproblem with respect to one of themodels and hence we do not need to readjust weights during minimization to accommodatethe different speed of convergence of the two models. Additionally, both inversion codesare altered only by adding the similarity term to the original inversion for a single typeof data. The drawback of the block descent is slower convergence (linear) towards theoptimal solution. This method is sometimes referred to as a Block Gauss-Seidel method.The algorithm is as follows:Following the Eq.(4.30), letting the Jacobians Je =∂u(σ)∂σ and Jf =∂s(m)∂m , and for simplenotation consider the measurement error matrices Σ−1e and Σ−1f to be unity, the GW andgeophysical descents can be summarized as53Algorithm 3 BCDMApproximately minimize Φ with respect to σApproximately minimize Φ with respect to mUpdate σ, s(m), if ∆m ≤ const. terminateDC descent step∇σΦ = JeTQe> (Qeu(σ)− de) + αeL>L(σ − σref) + β∇σΦS(s(m),σ)∇σ2Φ = Je>Qe>QeJe + αeL>L + β∇2ΦS(s(m),σ)∆σ = − (∇σ2Φ−1)∇σΦσk+1 = σk + µ∆σ (4.31)GW descent step∇mΦ = Jf>Qf> (Qfs(m)− df ) + αfL>L(m−mref) + βJf>∇mΦS(s(m),σ)∇m2Φ = Jf>Qf>QfJf + αfL>L + βJf>∇2ΦS(s(m),σ)Jf∆m = − (∇m2Φ−1)∇mΦmk+1 = mk + µ∆m (4.32)We still need to assign the weights α and β to balance the contribution of regularizationand similarity term. Weights α were set the same as for single data inversions and weightsβ were readjusted during the inversion process based on the data misfit Φd and ΦS :β = wΦdΦSwhere w ∈ [0, 1]. (4.33)The parameter µ, determining the size of the update for a given direction, is searched viasoft Armijo line search [75]. The convergence of BCDM where each block is not solvedexactly has been a subject of [83], more about the general convergence and conditions forthe sublinear rate can be found in [9].BCDM can be also applied if we know the petrophysical relationship, in such case thesimilarity term ΦS would be simply replaced by the petrophysical constraint. However, insuch case, the ADMM is expected to converge faster towards optimal solution than BCDM,and would be thus the preferred method.544.4.2 Joint Gauss-Newton method (JGN)In the context of this work, by Joint Gauss-Newton method we refer here to solving thecoupled problem jointly, updating both the groundwater and geophysical variable duringeach iteration. This means that all the terms of the objective function in Equation (4.30)are evaluated together. Hence compared to the BCDM discussed above, we keep bothdata misfit terms together in one objective function, i.e. no splitting to groundwater andgeophysical subproblems occurs. The advantage of this approach is overall lesser number ofiterations; however, we are now solving a larger system. Additionally, we need to determineall weights in equation (4.30) to balance the contribution of the different data misfits andthe regularization and coupling terms.At each Gauss-Newton step we then update the vector v that contains both conductivityσ and groundwater parameter m, v = (m,σ)>:(∆m∆σ)= −(∇σ2Φ ∇σ∇mΦ∇m∇σΦ ∇m2Φ)−1( ∇mΦ∇σΦ). (4.34)The sensitivities with respect to m and σ, which are in fact implemented as matrix timesvector routines, are likely to have a different computational time; however, in this setupcan not be solved in parallel. Therefore, despite the over all lower number of iterations thecomputational time for JGN might not be significantly better than for BCDM minimization.The off-diagonal terms of the Hessian matrix are coming from the structure similarityterm ΦS only, since we do not have an explicit petrophysical relationship and need to bederived analytically.∇σ∇mΦ = ∇σ(Jf>Qf> (Qfs(m)− df ) + αfL>L(m−mref) + βJf>∇mΦS(s(m),σ))= ∇σ(βJf>∇mΦS(s(m),σ))∇m∇σΦ = ∇m (β∇σΦS(s(m),σ)) (4.35)Schemes for both of the methods are depicted in Figure 4.1.55Figure 4.1: Comparison of the JGN and BCDM scheme. Left: The scheme for a singleiteration of the Joint Gauss - Newton minimization, the two data misfit terms ΦD,regularization terms ΦR and a similarity term ΦS are computed for the initial estimatesof m and σ, both models are then updated via a Gauss-Newton iteration (Eq.(4.34)).Right: BCDM minimization starts with the Gauss-Newton descent for the groundwatermodel m only (Eq.(4.32)) involving also the similarity term, the new GW model mk+1serves as an input for the geophysical Gauss-Newton descent (Eq.(4.31)), leading to anupdate for the geophysical model, σk+1, the next BCDM iteration starts with mk+1and σk+1.56Chapter 5Applications - synthetic examples5.1 Seawater intrusion - state estimation with ADMMThe first synthetic example follows a seawater intrusion scenario, where the variable densityflow brings more complexity to groundwater modeling, but due to the increased conductivityof seawater it is also an ideal target for most of the electromagnetic methods. Therefore,both groundwater and geophysical data are often collected to monitor seawater intrusions.Seawater intrusion (SWI) is a complex process that occurs naturally due to small dif-ferences in density between freshwater and saltwater. Depending on the hydrogeologicalsetting, seawater can enter the coastal aquifers through preferential flow pathways far intothe interior, or remain in a close proximity to the coast [4]. Increased groundwater ex-traction, reduced recharge into aquifers, and other human activities can cause the SWI topropagate further inland (Figure 5.1). To monitor the SWI and manage coastal aquifers,representative groundwater models need to be developed. Such models can then provideexplanations for saltwater origin in the area, and can be used to validate different pump-ing scenarios to manage the saltwater front propagation and future freshwater demands([29, 55, 61] or [87]).Numerous research studies have been published regarding mapping the seawater intru-sions progress by geophysical methods, especially large scale surveys are convenient for SWIdelineation ([19, 31, 72, 93] or [68]). Due to the complex character or governing equationsfor seawater intrusion process, the groundwater and geophysical model usually stay sepa-rate in most of the studies. Hence, in the coupled approaches the information is exchangedonly via some extra data or reference models.57Figure 5.1: Effect of increased pumping on the propagation of the SWI front. Source: A.D.Wenner, P.E. Jacobsen, L.K. Morgan - Understanding seawater intrusions [poster].2013For example, in the work of Herckenrath et al.[45], the geophysical estimates from 1DTDEM soundings served as extra observations for the GW model (after transforming viapetrophysical relationship - Archie’s law). In their coupled approach the GW model is usedto interpret the data and guide the geophysical inversion. A different way of transformingthe information from geophysical estimates was introduced in Beaujean et al.[8]. The ERTderived conductivities were transformed via Archie’s law to salt mass fraction estimates.These estimates were then filtered using a cumulative sensitivity based on squared Jacobiansand served as extra data for hydrological inversion next to groundwater salt mass fractiondata. The inversion was performed with PEST using a gradient based method, and thereforeit was not applicable to large scale problems.In the following synthetic example, we solve the inverse problem by minimizing bothtypes of data, and assume the knowledge of the petrophysical relationship, which providesthe constraint in the minimization. We use groundwater and geophysical data (well salinitydata end electrical potentials) collected only for a single instance in time, and we estimatethe unknown current and initial solute distribution based on this data. Since in this case weassume to know the petrophysical relationship, we can apply ADMM method for the jointminimization. Besides, we also show results for the direct substitution approach discussedat the beginning of section 4.3.The fact that we solve for the solute content-electrical conductivity only, with theassumption of at least an approximate knowledge of other GW parameters can be regardedas naive; however, the same framework can be established for GW parameters such ashydraulic permeability, external fluxes and other GW variables as long as appropriate58sensitivities are derived. We do not try here to estimate all of these parameters at once,since despite having two different sets of data, the inverse problem is essentially highlyill-posed, and most of these parameters vary over the entire domain. However, to test therobustness of the proposed minimization scheme we alter the “known” GW parametersused in the inversion such as hydraulic permeability and dispersion values, to see the effecton estimates.First, in section 5.1.1 we introduce the two cases of seawater intrusion scenarios, bothof them have a similar boundary conditions setup, but different permeability and thuslead to different SWI fronts. We also describe both groundwater and geophysical datasurvey designs. In section 5.1.2 we present the results of ADMM joint minimization andthe errors of the salinity estimates for a specific groundwater data survey; besides, we runmultiple inversions to see how the change of GW data locations affects the estimate errors.The results for altering the true GW parameters during the inversion are also presented insection 5.1.3. In the last section we discuss and summarize all findings from the quantitativeand qualitative point of views.5.1.1 Parametrization for the synthetic scenariosWe created two different model problems in 3D representing seawater intrusion; one witha homogeneous permeability field (Case 1) and one with heterogeneous permeability field(Case 2). The heterogeneous case is based on the field study at the Kidd2 site in the FraserRiver Delta in Richmond, BC [71], where the delta slope deposits confine the sandy deltaicdeposits and a seawater wedge enters from the river. In both cases the boundary condi-tions followed the Henry’s benchmark problem with hydrostatic pressure for the seawardboundary and freshwater inflow rate for the inland boundary. The actual parameter valuesare presented in the Table 5.6 including the external fluxes representing the pumping rates.For the initial “unknown” solute mass fraction distribution at time t0 we let the GWmodel run forward up to a certain time. During this simulation, a pumping well is placedin the southwest part of the area. Afterward, we altered the external fluxes, and a singlepumping well was placed in the north-east area while the freshwater inflow flux was de-creased. The GW simulation then ran from the initial state at t0 up to time t1 for 300days, with a time step 15 days. The “true” initial and final solute distributions for bothcases can be seen in Figure 5.3. The external fluxes are changed at time t0 so that the GW59GW model Heterogeneous case Homogeneous caseGrid 44x32x12 cells 44x32x12 cellsCell size 1 x 1 x 1 m 1 x 1 x 1 mPermeability kSilty sand kx = 2× 10−12m2, ky = 4.4× 10−11, kz = 2× 10−14m2 kx = 4.4× 10−11m2,Fine and medium sand kx = 4.4× 10−11m2, ky = 4.4× 10−11, kz = 4.4× 10−12m2 ky = 2.4× 10−11,Silty clay (Fig. 5.2) kx = 10−14m2, ky = 4.4× 10−11, kz = 10−17m2 kz = 1× 10−12m2Porosity φ 0.35 0.35Dispersion D 0.0032 m2/year 0.0032 m2/yearViscosity 0.001 0.001Freshwater density 1000 kg/m3 1000 kg/m3Saltwater density 1025 kg/m3 1025 kg/m3QGW , pumping rate up to t0 [x,y] = [8,14], 0.16 d−1 [x,y] = [8,14], 0.16 d−1QGW , pumping rate up to t1 [x,y] = [26,26], 0.13 d−1y [x,y] = [26,32], 0.13 d−1Geophysical modelGrid 50 x 38 x 12 cellsCell size 1 to 4 mArhie’s law m 1.7Background σ 0.0065 S/mTable 5.1: Parametrization for the seawater intrusion Cases 1 and 2.Figure 5.2: The geological layers for the heterogeneous case based on Kidd2 site in FraserRiver delta [71], schema for kz field.model, used in the coupled inversion, could not simulate the initial solute content from azero distribution. Moreover, changes to external fluxes (such as different pumping schemesor reduced discharge) are also likely to happen in real conditions.We collect both types of data only at time t1. For the GW sampling we have twotransects of wells (with spacing of 7 m) and 3 depth samples collected (depth = 4, 7 and11 m) in Case 1, and two depth samples (z = 5 and 9 m) for the heterogeneous Case2. The position of the transects was later altered to test different experimental designs,however, here we present in detail the case with west-east locations x = 16 m and x = 24 m,see in Figure 5.4. Gaussian random noise with standard deviation 0.05 was added to allmeasured solute fraction values. No hydraulic head or pressure data were used in the60Figure 5.3: Seawater intrusion scenario: Upper left and right: True models for initial andfinal solute distribution for Case 1; Bottom left and right: True models for initial andfinal solute distribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and 0.75 are displayed.coupled inversion.For the geophysical data, the simulated solute fraction at time t1 was converted throughArchie’s law into bulk electrical conductivity, and potentials were solved through the DCforward model described in Section 2.2.1. There are many different options for the electrodelayout and measurement scheme; the following one was chosen based on the sensitivities ofmeasured data, while trying to maximize the depth resolution for data collected only at thesurface. The electrode layout corresponds to a regular grid with spacing 3m in x directionand 4m in y direction, giving in total 72 electrodes. A positive electrode was fixed closeto the seaward boundary (west) and the negative charge was moving along the x profile,61towards east. For each source pair (72 in total) potential differences were measured on allreceivers, where one of the receiver couples was always fixed and placed in the north westcorner (see Figure 5.4). 3% Gaussian random noise was added to the measured potentials.Figure 5.4: Data collection. Left: Experimental design for the DC survey: Dark blue pointsrepresent the electrodes on the surface placed on a regular grid, saltwater is coming infrom the west boundary. Right: GW wells locations, 8 wells in total are placed alongtwo transects at x = 16 m and x = 24 m distance, later labeled as b and c, otherexperimental designs for GW data are plotted in Figure 5.10.5.1.2 Joint inversion with ADMMThe ADMM minimization starts with the GW model descent and continues as long asthe constraint residual rk decreases (or up to 5 runs of GW and geophysical descent).The residual rk represents the difference between bulk electrical conductivity and electricalconductivity derived from the GW model via Archie’s law. Since this is a synthetic examplewe can record the actual initial and final errors next to the data misfits for both thegroundwater and geophysical data during the minimization. By actual error we mean thenorm (ωk) = ‖ωk − ωtrue‖. Due to the ADMM approach we do not need to weight thetwo different data misfits, however, weights still need to be assigned for the regularizationterm α and the so called penalty term, ρ (see the augmented Lagrangian formulation in(4.24)). The choice of the regularization parameter α has already been largely discussedin the literature (see [41, 94] and reference within), and can be determined either basedon initial values for φD, φS and φR, or by a trial and error procedure. The choice of thepenalty parameter ρ was discussed in the section 4.3.1. A particular set of weights wereapplied to obtain all results presented here in; αe = 10−3, αf = 5 × 10−3 and ρ = 0.4.62The number of Gauss-Newton iterations (within each ADMM descent) was 3-4 for the GWblock descent, and between 6 to 10 for the geophysical block descent.In the Figure 5.5, the actual errors scaled against the error of initial estimates areplotted together with the residual rk as they decrease during the minimization. Initialestimates are based on the forward simulation starting with the GW reference model. Bothrk and k decrease during the ADMM minimization. The estimates of ω0 can be seen inthe contour profiles in Figure 5.6, resp. Figure 5.7 for x = 10 and x = 20 m or their 3Dplots in Figure 5.8.Case 1 Case 2Figure 5.5: The error and residual decrease during the ADMM descents: Green trianglesrepresent the scaled error for the final solute fraction ωf , orange stars correspond toupdated rk values, where rk = σf (k) − p(ωf (k)). The GW wells in this case wereplaced along x = 16 m and x = 24 mFor comparison, we also solved the inverse problem with a simpler coupled approachand by a Gauss-Newton method with a direct substitution (equation (4.23)). In the simplecoupled approach both models can run independently. We first solve the inverse problemfor GW data only and then apply Archie’s law to transform the estimate to electricalconductivity at t1. This estimate then constrains the geophysical inversion as a referenceand initial model, and as such it is computationally easier to implement with no extracoupling terms in any of the objective functions. The actual errors of solute mass fractionat t0 and t1 and final data misfits for the ADMM and the coupled approach are presentedin Table 5.2. In Table 5.2 we recorded for the coupled approach the lower error from GWor geophysical inversion. With the ADMM approach we decreased the error for the finalsolute fraction ωf by roughly 50% compared to the coupled approach, and also decreased63Figure 5.6: ADMM estimates for Case 1: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full contourlines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5 (green) andω = 0.75 (red).Figure 5.7: ADMM estimates for Case 2: Contour profiles at x = 10 and x = 20. Thedashed lines are the estimates from the joint inversion (ADMM), and the full contourlines are for true locations corresponding to ω = 0.25 (blue), ω = 0.5 (green) and ω= 0.75(red).the error for the initial estimate ω0.Direct substitutionThe Gauss-Newton minimization with a direct substitution for σ follows the objective func-tion in equation (4.23). We run the Gauss-Newton minimization applying three differentstrategies of adjusting the weights and summarized the results for Case 1 in Table 5.3. In64Figure 5.8: ADMM estimates. Upper left and right: Estimates for initial and final solutedistribution for Case 1. Bottom left and right: Estimates for initial and final solutedistribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and 0.75 are displayed. The truemodels are plotted in Figure 5.3Case 1 (ωf ) (ω0) φGW (ρω) φDC(u(σ))Initially 14.2 22.6 14.4 1046Coupled 6.28 13.56 4.6 1.6ADMM 3.1 10 1.1 1.4Case 2Initially 9.3 16.8 7.4 1940Coupled 6.6 13.8 3.4 1.0ADMM 2.8 11.5 0.8 0.9Table 5.2: Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 for the twodifferent cases.65the first option we initially fixed both δe and δf to make the initial data misfit magnitudesequal. In the second option we fixed the geophysical data misfit weight δe, and adjust thegroundwater data misfit weight δf at each iteration as δf =φD,eφD,f. In the third option wewere adjusting the geophysical weight δe instead, at each iteration as δe =φD,fφD,e. Since theestimates by substitution approach were worse than computationally much simpler coupledapproach described above, we left the substitution approach out of further comparison.Case 1 (ωf ) (ω0) φGW φDCInitially 14.2 22.6 974 1046Fixed δe, δf 10.8 27.4 7.97 43.74Adjusted δf 8.49 21.3 1.78 144.24Adjusted δe 7.07 20.6 3.2 71.12Table 5.3: Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 for Gauss-Newton method with the substitution of σ; 1. The weights were fixed throughout theminimization based on the initial data misfits; 2. δf was adjusted as δf =φD,eφD,fat eachiteration; 3. δe was adjusted as δe =φD,fφD,eat each iteration. GN inversion run up to 10iterations.Different GW experimental designsAdditionally, we tested the ADMM and the coupled approach for different locations ofGW wells without changing the DC survey design. The scaled errors are plotted for allsimulations in Figure 5.9. The different transects of wells are plotted in Figure 5.10. Wedid not use the same combinations of transects for Case 1 and 2, as the final SWI frontreached further in the Case 2 compared to Case 1. In most of simulations we managed tosubstantially reduced the error for final solute mass fraction ωf (about 25% of the initialerror for homogeneous Case 1 and 35% for the heterogeneous Case 2). The errors for initialsolute estimate were fairly consistent and satisfactory for the homogeneous Case 1 (50%reduction), less so for the Case 2 (60 to 80%).5.1.3 Joint inversion by ADMM with inexact GW parametersIn order to test our method for the case, where the reservoir parameters are not knownexactly or only approximately, we solve the problem for an inaccurate permeability field anddispersion. In the first test, we used the original GW model parameters, but decreased thehomogeneous permeability field and dispersion to 70% when running the ADMM or coupled66Case 1 Case 2Figure 5.9: The relative error decrease for five different GW sampling designs; the dottedline is for the ω0 relative error, the full line for the ωf relative error decrease. In allcases the error decrease slows down with further iterations. The plotted results arebased on different transects of wells plotted in Figure 5.10.Figure 5.10: The transects along x and y axis in plan view. For the homogeneous Case1the sampling depths were z = [3, 7, 11], and for the heterogeneous Case 2 z = [5, 9].each design had two vertical or horizontal transects.67approach. The ADMM joint approach converged, but the actual errors were higher thenwhen the correct GW parameters are used. In Figure 5.11 you can see the error evolutionsfor both the ADMM with correct and incorrect GW parameters, also the error (ωf ) fromthe coupled approach (a single value), Table 5.4 provides the summary of errors for bothmethods.In the second test, we used the homogeneous Case 1 and altered the permeability field byadding a 3D random Gaussian field to the original one (see in Figure 5.12). The addition ofthe Gaussian random field thus changed the original anisotropic homogenous permeabilityfield to a heterogenous field, which leads to a differing solute distribution. The ratio ofchange in the observed solute fraction data due to different permeability compared tooriginal data was 18%. GW data based on this simulation where used in ADMM inversion,but leaving the permeability field homogeneous as in the previous calculations for Case 1.For comparison we again ran the coupled approach with the same input as used in ADMM.The scaled errors and residual rk are plotted in Figure 5.13 for both correct and incorrectGW parameters.Figure 5.11: Altering the GW parameters: Errors decrease for ωf estimate and residual rkwhen correct and altered GW parameters (D, k) are used in ADMM. The errors forωf for the coupled approach are plotted as well with a dashed line, it is just a singlevalue.68Figure 5.12: Altering the GW parameters: The true permeability field in x direction kx inthe second test. For solving the inverse problem a fixed value kx = 4.4−11 was usedinstead, while the true field is shown in this figure. The true ky and kz were alteredsimilarly when simulating GW data.Figure 5.13: Altering the GW parameters: Errors decrease for ωf estimate and residual rkwhen correct and altered GW permeability is used in ADMM. The errors for ωf forthe coupled approach are plotted with a dashed line.69Test 1 (ωf ) (ω0) φGW (ρω) φDC(u(σ))Initially 12.2 20.7 63 1150Coupled 7.35 16.2 3.64 2.18ADMM 3.85 10.4 2.04 1.57Test 2Initially 17 22.6 17.2 1140Coupled 6.92 13.1 4.0 1.95ADMM 4.32 9.26 1.7 1.48Table 5.4: Errors, (ωk) = ‖ωk − ωtrue‖, for the solute content at time t0 and t1 whendifferent from true GW parameters are used in the ADMM inversion or coupled ap-proach. Test 1 - change in permeability field. Test 2 - 70% reduction in permeabilityand dispersion values.705.1.4 SummaryWe can conclude that the ADMM approach proved to be advantageous compared to thesimple coupled approach, direct substitution of the constraint or the separate groundwaterand geophysical inversions. In all tested examples, the joint inversion with ADMM achieveda lower error for both the initial and final solute fraction distributions compared to thesimple coupled approach.More particularly:• Based on the different examples presented herein the error for ωf estimate by ADMMwas roughly 50% of the error by a coupled approach in homogeneous case, and 60%for the heterogeneous case. For the initial solute fraction ω0 the improvement byADMM was 60% and 70% respectively.• The estimates for the final solute distribution were generally better than for theinitial solute mass fraction due to the fact that the DC data were collected at thefinal time, and also the coupling constraint between ωf , solute fraction, and σf ,electrical conductivity was enforced for the final time t1.• In all cases the ADMM converged to minimum, though it shows some of its typicalaspects: a relatively quick drop during the first few iterations and a slow decrease to-wards the end. However, for the case of groundwater modeling applications, the initialdecrease in error might be sufficient. Our synthetic study confirmed the theoreticalresults about the convergence of this method.Of course, the ADMM comes at a higher computation cost compared to the simplecoupled approach, which is due to the repeated computation of the GW and geophysicaldescent, for four or five times. For example, using the simple coupled approach for the Case1, we ran the forward GW model 10 times, while the ADMM required 57 forward modelruns. Similarly, the coupled approach required 7 DC forward model runs as opposed to 32runs by ADMM. In the Table 5.5 we summarized the computation cost in terms of forwardgeophysical and groundwater model runs and the actual running time; we also included thetime scale of the Gauss-Newton method with a direct substitution.Note that, when the relative weighting and regularization parameters are unknown, theADMM can be easily faster than the Gauss-Newton optimization (with substitution forω0), which also needs to be solved a number of times due to testing different regulariza-tion parameters, however, unlike ADMM, where each iteration consists of solving a single71physical model, the fully coupled Gauss-Newton method requires solving both problems foreach weighting and regularization parameters.Method GW f.m. runs DC f.m. runs time [min]Substitution (8 iterations) 20 18 59Coupled (4 and 7 iterations) 10 7 27ADMM (5 descents x 3 and 6 iterations) 57 32 113Table 5.5: Computational cost overview for the forward model (f.m.) runs and actual runtime of the inversion; all results displayed were executed on a desktop with Intel(R)Core(TM)i7-2860 QM processor and 16GB RAM. Parameters for line search procedurewere the same for all methods, amount of iteration counts is displayed in the bracketsnext to each method.For the first set of simulations with Case 1 and Case 2 we assumed the GW modelparameters were known, excluding the initial solute content. This can be regarded as anoverly simplifying approach, but is justified for testing the feasibility of the joint inversionstrategy. To demonstrate the robustness of the ADMM method with respect to variationsin GW parameter values, we altered the GW model parameters in the inversion process.As expected, this led to estimates with higher error compared to solving the problem withthe correct parameters. Nevertheless, the ADMM converged to estimates with lower errorthan the simple coupled approach (see in Figure 5.11), as it could be partially “corrected”by information from geophysical data. We are aware that further increasing of the error inthe GW model parametrization would also lead to worse estimates, but in that case anycoupled approach will not be able to provide a more accurate estimates.725.2 Solute tracer test - Parameter and state estimationwithout an explicit petrophysical relationshipThe solute tracer test provides a relatively simple case, on which we can test using implicitstructure similarity measures. Abundant research has been published on the hydrogeophys-ical inversion for solute tracer tests, for both field and synthetic studies; however, only fewstudies have been published which did not use an explicit petrophysical link.For example in [13] authors implemented a level set method in the objective functionand assumed the parameter of interest is mainly affected by geological facies. The approachwas tested with seismic and hydrological tomography data on a synthetic example. Dueto the level set method there was no need for a direct link between seismic velocity andhydraulic conductivity; however, in this application the structure similarity was enforcedby the a priori assumption, i.e. the knowledge of geological facies location, and could beapplied only due to the fact that both geophysical and hydrological data were directlydependent on the same geological property.A different way to alleviate the need of petrophysical relationship is to assume a strongcorrelation between the variables of interest, in the study of [54] between the solute concen-tration and soil bulk conductivity. In particular in their synthetic study for a solute tracertest, the authors used ERT time lapse sampling together with groundwater head and fluidconductivity data to determine the hydraulic conductivity. In their approach they maxi-mized the correlation between transient changes in the fluid and soil bulk conductivity andthus avoid solving the ERT inversion; the input for the ERT forward model were fluid con-ductivity changes alone. The correlation “coupling” term was then added when solving thestandard groundwater inverse problem to estimate hydraulic conductivity. The hydraulichead and fluid conductivity sensitivities were solved by perturbation method in paralleland hydraulic conductivity was parametrized by Pilot points method with a linear interpo-lation; a high performance cluster was necessary to run the inverse problem. The presentedexample showed that even indirect geophysical information can improve the estimates forhydraulic conductivity.To our knowledge, the first work applying a structural constraint in hydrogeophysicalinverse problem is the study of Lochbuhler et al.[65]. The authors combined two datasets,ground-penetrating radar data with hydraulic tomography or tracer mean arrival times, toestimate hydraulic conductivity. A cross product of the spatial gradients of the geophysicaland groundwater models was implemented as a structural constraint, creating a penalty73in the coupled objective function. Their study tested the structural coupled approach on2D synthetic case as well as field data from the Widen gravel aquifer in Switzerland. Therecovered models did indeed displayed similar structure, though not following a simple linearpetrophysical relationship, and the hydraulic conductivity estimates were in accordancewith previous field studies.In our synthetic case, we implemented a full 3D inversion and jointly inverted the geo-physical and hydrological data that are not directly dependent on the property of interest.We used the Cross Gradient field product as a structure similarity coupling term in theobjective function in Eq. 4.30 and solved the inverse problem for either initial solute stateor hydraulic conductivity field. Two computationally different methods, Block CoordinateDescent (section 4.4.1) and Joint Gauss Newton (section 4.4.2) were tested and comparedwith the results of separate geophysical and groundwater inversions.5.2.1 ParametrizationWe set up two different concentration plumes, first one is initially a standard Gaussianplume (Example 1), and the second one has a more complex shape (Example 2). Forthe first two examples we estimate the initial conditions for solute concentration, c0, inthe last case (Example 3) we used the same initial concentration as in Example 1 butestimate the hydraulic conductivity K instead. For each case, the boundary conditions areset either no flow or fixed head boundary, such that the natural hydraulic head gradientis along the x coordinate. Injection and pumping wells are placed in the domain to createflow pathways also along the y and z axis. The solute plume moved from time t0 totime t1, and the solute content at t1 was then transformed via Archie’s law to create anelectrical conductivity model, and as such used in the Example 1. To imitate more realisticcases, in the subsequent examples the electrical conductivity model was altered by addinga heterogeneous spatially correlated field of more conductive material (see in Figure 5.14)than background. The actual parameters for both the groundwater and geophysical modelscan be found in Table 5.6.Both groundwater and geophysical data were collected at a single instance of time, t1,though multiple time collection would not change computationally the algorithms presentedin chapter 4. The groundwater measurements were collected from a few nested wells, toobtain solute concentration data. For the Example 3 we also recorded the hydraulic headdata.74Table 5.6: Parametrization for the solute tracer Examples 1, 2 and 3.GW modelGrid 42x32x16 cellsCell size 1 x 1 x 1 mHydraulic conductivity KExample 1,2 kx = 1m/d, ky = 0.6m/d, kz = 0.4m/dExample 3 kx,y,z , log(µ) = −1.8m/d, σ = 1, lx,y,z = (60, 40, 10)Dirichlet BC along x h(x = 0) = 21.6 m, h(x = 42) = 20 mPorosity φ 0.35Dispersion D 0.02 m2/dQGW , injection rate [x,y] = [8,8], 4 d−1QGW , pumping rate [x,y] = [34,18], 3.3 d−1Geophysical modelGrid 50 x 38 x 16 cellsCell size 1 to 4 mArhie’s law m 1.7Background σ 0.0065 S/mFigure 5.14: Solute tracer test. Left: A scheme for the GW survey design with twotransect of wells, applied in Example 1 and 2. Right: A scheme for DC surveywith receivers on surface only, in the background is the electrical conductivitymodel applied in Example 2 which spatially correlated more conductive fieldthan background.The placing of the wells has a huge impact on the inverse problem solution and itis generally hard to design without the actual knowledge of the solute plume location;therefore, we ran all the examples for different groundwater experimental designs, withwells randomly chosen from a regular grid of 4 x 5 sampling locations. All different GWsurveys wells locations are plotted in Figure 5.15.For the geophysical survey, the potentials were measured using surface electrodes only.A regular grid spacing was 4 m and 3m along the x and y coordinate. The source electrodecouple was moved along x (north - south) coordinate and for each source the potentials75Figure 5.15: Different GW surveys for solute tracer test. Ex 1,2 - des1: Two transects ofwells at x = 14 and 18m applied in Example 1 and 2. Ex 1,2 - grid: 4 x 5 grid ofwells from which 8 locations was chosen randomly for each run. Ex3-des1: 3 x 3grid of wells locations for results presented in Table 5.11. Ex3-grid: Regular grid ofwell from which 10 wells were randomly chosen, corresponding results summarizedin Table 5.12GW model Geophysical modelWeights α 1e− 3 1e− 4β0 27 27β 0.5 φD/φS 0.5 φD/φSIterations separate inversion 4 6BCDM, GW / DC max 3 4JGN, 7 7Conjugate gradient 8 8Table 5.7: Minimization parameters for the inverse problemwere measured by all receivers with one electrode being fixed and located in the south-eastcorner. The design follows the one already used in the seawater intrusion example.For all tested scenarios we first ran separate inversion for groundwater and geophysicaldata alone, and afterward, we apply two computational approaches for the joint inversion,Gauss-Newton method (JGN) to update both the groundwater and geophysical models atonce, and BCDM with up to four groundwater or geophysical descents. Each inversionwas stopped when the data misfit dropped bellow expected error or when the maximum ofiterations was reached. The iterations and weights setup are summarized in Table 5.7.To compare the different methods for parameter or state estimation, we computedthe norm of errors (since the true parameters are known) for the initial and final solute76estimates, soil bulk electrical conductivity and in the last example also for the hydraulichead and conductivity estimates. We also evaluated the Cross Gradient similarity termbefore and after the structure coupled minimization. The visual comparison is somewhatlimited as all the examples are in 3D, however, vertical and horizontal slices are presented.5.2.2 Example 1 and 2: Estimation for initial state c0Figure 5.16: Solute tracer test. The true models for the initial solute concentrationin Example 1 (left) and Example 2 (right), isosurfaces displayed at c0 =0.2, 0.4, 0.6 and 0.8.First, we applied the inversion framework to a simple case of a solute tracer test wherethe initial solute plume has a Gaussian shape (see in Figure 5.16). The true solute content cnis related to electrical conductivity σ via Archie’s law. Solute concentration and geophysicalpotentials were collected at time tn. In the Table 5.8 we summarized the results of separateinversions and joint inversion with block coordinate method (BCDM) and direct Gauss-Newton method (JGN) for a particular placing of the GW wells (marked in Figure 5.15 asEx1,2-des1). We recorded the final error estimates, data misfits and a drop in the similaritypenalty term, cross gradient product φCGP . These errors are visualized in the bar plot inFigure 5.17, the separate (GW data only) inversion decreased the initial error of cn to 60%while BCDM decreased the error to 45%. The joint inversion solved with JGN had worseperformance, with 55% error decrease only. The errors for the estimate of initial solutediffered less among the applied methods, separate inversion had 72%, BCDM 65% andJGN 70%. Although all results were computed in 3D, in the Figure 5.18 we show the 2Dcontour plots of solute estimates for cn along y = 16 and depth slices at z = 8, the initial77solute estimates are plotted for z = 8.Method (cn) (c0) φD(cn) (σ) φ(σ) ΦCGP (cˆn, σ) ΦCGP (cn, σˆ)GWI 5.20 7.06 2.71 - - 24.9 -DCI - - - 4.93 2.16 - 34.0JHI (GNm) with φCGP 4.81 6.88 8.12 4.18 7.95 19.03 2.88BCDM with φCGP 3.94 6.35 1.15 4.32 16.9 17.05 1.46Table 5.8: Example 1: Summary of the errors for the solute distribution and electricalconductivity estimates, the data misfits and the cross-gradient coupling terms (betweenthe estimate and the true field) by either separate inversion or joint inversion withcross-gradient product.Figure 5.17: Example 1 and Example 2: Bar plots for the error estimates by differentmethods, the actual values are in Table 5.8, Table 5.10 respectively.Error cˆn cˆ0 σˆExample 1 mean (JGN/GWI/DC) 0.88 1.25 0.98Example 1 mean (BCDM/GWI/DC) 0.84 0.91 1.00Example 2 mean (JGN/GWI) 0.91 0.99 1.1Example 2 mean (BCDM/GWI) 0.78 0.84 1.11Table 5.9: Example 1 and 2: Summary of the relative errors (joint inversion vs GW (or DC)data only) for the solute distribution and electrical conductivity estimates based on theresults from 10 different GW survey designs.For the Example 2 we altered the true initial solute distribution by adding a spatiallycorrelated field of higher conductivity (see in Figure 5.16). Unlike in the previous example,the corresponding soil bulk conductivity is therefore not linearly dependent on solute con-centration model. The summary of the estimate errors by different methods can be foundin Table 5.10 and the slices of the actual estimates are plotted in Figure 5.19 for a particular78Figure 5.18: Example 1: The true models vs the estimates by GW data inversion andstructure-coupled inversion, contours at z = 8 and y = 18 for final and initial solutecn, c0 are plotted.well placing. Similarly as in Example 1, BCDM had lower errors than JGN, decreasing theinitial error to 43% for the solute content cn, and to 54% for c0. The separate inversion withGW data only decreased the error to 58% for cn, respectively 70% for c0. The estimates forelectrical conductivity were comparable regardless the method, only BCDM with Archie’slaw lead to significant estimate error reduction.Afterward, we ran both Example 1 and 2 for different GW survey designs; each timeeight wells locations were randomly picked from a regular grid covering the area of plumemovement. In Table 5.9 we summarized the results of these experiments, the actual errorsfor each run of Example 2 setting are plotted in Figure 5.20. On average, BCDM reducesthe error of the final solute content by 16 % in Example 1, and by 22 % in Example 279compared to GW data inversion alone. For the initial solute the drop was 9% and 21%respectively. JGN also reduced the final cn estimate errors, but has not improved the initialsolute estimate compared to GW data inversion.Method (cn) (c0) φD(cn) (σ) φ(σ) ΦCGP (cˆn, σ) ΦCGP (cn, σˆ)GWI 7.03 9.43 0.9 - - 108.3DCI - - - 6.85 6.4 - 32.6JHI (GNm) with φCGR 6.29 8.53 6.8 7.04 7.1 73.9 9.4BCDM with φCGR 5.18 7.29 0.9 6.93 7.3 32.6 6.9Table 5.10: Example 2: Summary of the errors for the solute distribution and electrical con-ductivity estimates, the data misfits and the cross-gradient coupling terms (betweenthe estimate and the true field) by either separate inversion or joint inversion withcross-gradient product.5.2.3 Example 3: Estimation for hydraulic conductivity KThe 3D variable hydraulic conductivity K is usually parametrized so that its dimensionis reduced. Multiple approaches, which can be generally cast as a linear combination ofsome M basis vectors (where M is usually much smaller than the mesh size) can representthe hydraulic conductivity K [76]. The simplest case is a zonation technique, where thesebasis vectors are constant over a part of the domain and zero everywhere else. Anothercommon option is a pilot point technique; the hydraulic conductivity values are fixed atsome points and the rest of the domain is determined by krigging interpolation. The choiceof the locations of these points, or adding additional ones during the minimization leads todifferent approaches, examples can be found in [82, 86].For Example 3 we assumed that the hydraulic conductivity field is log-normally dis-tributed with some average value µ and covariance C.ln(K) ≈ N(µ,C)The covariance C is usually derived assuming Gaussian or exponential correlation lengthsin the x, y and z direction and some known variance (sometimes referred to as a semi-variogram). The semivariogram values are mostly obtained from geology borehole logs orother available data. Such geostatistical conceptualization is widely used, for example in[51] or [79]. Given this geostatistical interpretation, it is possible to apply Karhunen-Loeve80Figure 5.19: Example 2: The true models vs the estimates by GW data inversion andstructure-coupled inversion, contours at z = 8 and y = 18 for final and initial solutecn, c0 are plotted.decomposition and write the covariance for the discrete hydraulic conductivity field asC =N∑j=1λjvjv>j , (5.1)where C is an N × N discretization of the covariance operator, C, the vectors, vj are itseigenvectors, that is, v>i vj = 0 if i 6= j and 1 otherwise and λj are the (positive) eigenvalues.When solving the inverse problem, we do not need to work with a full decompositionof the covariance matrix and rather replace it with the first M basis vectors,ln(K) = µ+M∑j=1wjmj = µ+ Vm (5.2)81Figure 5.20: Example 2: The estimates for 10 different experimental designs, altering loca-tion of 8 wells from Ex 1,2 - grid design.where wj =√λjvj and V = [w1, . . . ,wM]. M should be large enough to capture theheterogeneity of the conductivity field K. This significantly reduces the size of the inverseproblem to estimating the weights of these M basis functions, stored in the vector m;however, the drawback is that we limit any solution to a Gaussian random field given bythe covariance matrix C.The unknown conductivity field K was thus parametrized asK = exp(µ+ Vm), (5.3)and the covariance matrix in this example was computed based on Gaussian exponentialmodel:Ci,j = −12σ2 exp((xi − xj)2υ), (5.4)where υ corresponds to correlation lengths in the x, y or z direction, σ is variance and(xi − xj) is the distance between the two points. The values of vector υ and σ2 are laterused to numerically solve the Karhunen - Loeve expansion of the covariance matrix C up82to the first M vectors.The number M can be estimated from the plot of eigenvalues, and verified by simplysolving a linear least square problem with the actual data for ln(K):mˆ = (W>W)−1(W>ln(K)),where W is the matrix containing the basis vectors from the covariance matrix decompo-sition and a unit vector e representing the mean µ; W = [e V]. The “reduced”hydraulicconductivity K(mˆ) can be then compared against the true hydraulic conductivity K. Inour case M ranged in the interval (30,60) to sufficiently describe the heterogeneity of thehydraulic conductivity, however this number is strongly dependent on the vector of correla-tion lengths υ. The decomposition of the matrix C was solved numerically using Choleskyfactorization with Matlab built-in routines. Lancozs algorithm would be an alternative forlarge scale problems.The setup of Example 3 is the same as in previous cases except that now we have aheterogeneous hydraulic conductivity K and hydraulic head data were collected next tothe solute concentration data, providing indirect information about hydraulic conductiv-ity. First, we placed nine wells in regular spacing (see in Figure 5.15) and the inversionparameters were initially set as zero for the m weights, and initial mean value µ0 = −1and M = 35. We again recorded the estimate errors for final solute content and electricalconductivity, but this time also for hydraulic head, µ and ln(K) in the Table 5.11.Afterward, we run the inversions for different groundwater sampling designs, choosingrandomly 10 wells from the well points grid (5.15). The summary of these results can befound in Table 5.12 and the actual errors are plotted in the Figure 5.23. For all results thevector m was set initially as zero vector, and the initial mean value as µ0 = −1.In summary, the closer was the initial µ estimate to the true value, the better estimates(lower actual error) of solute distribution were obtained. The errors of the estimates ofhydraulic conductivity and also final solute cn were lower for BCDM and JGN comparedto GW data inversion alone, however not in all cases BCDM outperformed uncoupledinversion and in some cases the final estimates were fairly comparable. On average JGNdid not manage to minimize the hydraulic head data misfit, and BCDM had thereforeoverall better performance.83Method (cn) (h) (σ) (K) µˆ C(cˆn, σ) C(cn, σˆ)GWI 5.18 36.8 - 78.9 -2.26 50.41 -DCI - - 3.37 - - - 6.32JHI (GNm) with φCGR 4.41 38.7 4.56 68.6 -1.54 26.99 4.32BCDM with φCGR 3.52 20.1 3.64 63.6 -1.93 35.23 0.85Table 5.11: Example 3: Estimation of hydraulic conductivity with KL expansionparametrization: Summary of the errors for the solute distribution, electrical conduc-tivity and hydraulic conductivity estimates (and also the estimate for µ = mean(K),the true value was mu = −1.8) by either separate inversion or joint inversion withcross-gradient product. A final cross-gradient product between cn and σ is displayedtoo.Figure 5.21: Example 3 : Bar plot for the errors by different methods, the actual values arein Table 5.11Error cˆn hˆ ln(Kˆ)Example 3 mean (JGN/GWI) 0.88 2.15 0.86Example 3 mean (BCDM/GWI) 0.77 0.92 0.8Table 5.12: Example 3: Summary of the relative errors (joint inversion vs GW data only) forthe solute distribution at cn, hydraulic head h and logarithm of errors for conductivityK based the results from 10 different GW survey designs.84Figure 5.22: Example 3: The true models vs GW data and structure-coupled inversionestimates for the final solute concentration, contours at z = 8 and y = 16 are plottedand hydraulic head at y = 24.85Figure 5.23: Example 3: The estimates for 10 experimental designs by different methods,altering the location of 10 wells in Ex3-grid design.86Figure 5.24: Example 3: The true field and estimates of hydraulic conductivity K by GWdata only and JGN and BCDM.875.2.4 SummaryBoth methods, separate data inversion and joint approaches, were able to minimize thedata misfits and provide physically realistic estimates for both hydrological and geophysicalmodels. The coupled framework (joint inversion) reached better estimates for both initialsolute content and hydraulic conductivity in most of the scenarios presented above andthus gave a promising improvement without a need for an explicit petrophysical link.In particular, from evaluating the numerical results and comparing the actual estimateswe can conclude the following points:State estimation (Examples 1,2)• BCDM had the lowest errors for both the initial and final solute estimate;• JGN had lower errors for the final solute estimate in all cases compared to GW datainversion alone, however, for most experimental designs the initial solute estimatewas worse than in GW data inversion alone;• Soil bulk conductivity estimates were not improved by the joint inversion (both JGNand BCDM);• The more complex Example 2 (no linear relationship between soil conductivity andsolute content) lead to improvements for both BCDM and JGN methods comparedto GW data inversion alone.Parameter estimation (Example 3)• JGN and BCDM provided comparable results in terms of (K) and (c);• The initial estimate for µ does influence the final errors;• In the hydraulic head estimation JGN performed worse than both GW data inversionand BCDM;• The main features of the heterogeneous hydraulic conductivity K were captured byboth coupled and uncoupled inversion framework, however GW data inversion con-verged towards higher contrast in K.88• To our knowledge, no structural inversion have been done using a parametric ap-proach. Here, we explore the use of parametric approach together with a structuralcoupling.The GW inversion is mainly limited by the scarce groundwater data, and as such it couldnot provide better results without implementing some smooth regularization, however, inthat case, “a-priori” information needs to be introduced, which can not be always justified.Adding a coupling term in a form of structure similarity measure better constrains the ill-posed GW inverse problem and complements the GW data. This observation also partiallyexplains why in Example 2 BCDM provided on average better estimates than in Example1. Example 1 was set as a simple Gaussian plume, and thus standard L2 regularizationwithin GW inversion alone gave satisfying results, however, if more complexity is added(Example 2), then, the geophysical contribution within the CGP term contained sufficientinformation to constrain the estimates.Both JGN and BCDM were easier to implement for the state estimation of the solutecontent c0, due to easier derivation of the data sensitivities. The coupled approach for stateestimation had better results when compared with the coupled approach for parameterestimation. This stems from the fact that the similarity measure term is coupling thesolute content c with the electrical conductivity σ. Other GW parameters such as hydraulicconductivity or, say, boundary fluxes are related to solute content c via the GW flow modelonly and this dependency has then a major effect on the success of parameter estimation.895.3 Joint inversion with electromagnetic loop surveyThe presented results follow the seawater intrusion example, except that now we haveTDEM data instead of DC resistivity data. The details of this example were alreadydescribed in section 5.1, and here we want to test the robustness of the ADMM approach ifthe geophysical data are altered. The groundwater model boundary conditions and forwardsimulations are identical, and again we estimate the initial and final solute distribution forthe two cases with the homogeneous and heterogeneous permeability field, Case 1 and Case2 (see Table 5.6).To model geophysical electromagnetic data we placed a transmitter loop at the centerof the area on surface. The current was initially shut down to produce changes in magneticfield, which were then recorded by receiver coils placed in a uniform grid on surface (seeFigure 5.25 for the experimental design scheme). 3% Gaussian random noise was added tothe measured data.Figure 5.25: The experimental design for EM survey: Loop (transmitter) is placed in thecenter on surface, the receiver coils are placed on surface in an uniform grid.For the GW experimental survey we again used the two transects at x = 16 m andx = 24 m, with four wells each, to collect solute mass fraction data from three depthsamples.905.3.1 Joint inversion with ADMMThe ADMM starts with a groundwater descent, and continues as long as the constraintresidual rk decreases. The weights were set as αe = 10−9, αf = 5 × 10−3 and ρ = 0.4.The summary of results for the case with well transects along x = 16 m and x = 24 m isin Table 5.13, Figure 5.26 shows the scaled errors (k/0) and the residuals rk during eachdescent of ADMM method for this particular survey. Similarly as the results in section 5.1showed for DC resistivity data, also when coupling groundwater with electromagnetic datathe ADMM had the lowest errors for both initial and final solute mass fraction estimates.Case 1 (ωf ) (ω0) φGW (ρω) φDC(u(σ))Initially 14.2 22.6 14.4 1.2 e5Coupled 6.1 15.3 5.1 70.8ADMM 2.6 13 2 69Case 2Initially 9.3 16.8 7.4 7.1 e5Coupled 5.2 14.1 2.9 75.7ADMM 2.0 10.4 1.8 71.3Table 5.13: Errors for the initial and final solute content at time t0 and t1 for the twodifferent Cases 1 and 2.Case 1 Case 2Figure 5.26: The error and residual decrease during the ADMM descents: Green trianglesrepresent the scaled error for the final solute fraction ωf , orange stars correspond toupdated rk values, where rk = σf (k) − p(ωf (k)). The GW wells in this case wereplaced along x = 16 m and x = 24 mThe profiles of initial and final solute mass fraction estimates versus the true modelsare plotted for Case1 and Case 2 in Figure 5.27, Figure 5.28 respectively.91Figure 5.27: ADMM estimates for Case 1: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full contourlines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5 (green) andω = 0.75 (red).Figure 5.28: ADMM estimates for Case 2: Contour profiles at x = 10 and x = 20. Thedashed lines are estimates from the joint inversion (ADMM), and the full contourlines are the actual locations corresponding to ω = 0.25 (blue), ω = 0.5 (green) andω = 0.75 (red).Different GW experimental designsAfterward, we run the joint inversion for different locations of GW data. We used the samecombinations of transects as before, plotted in Figure 5.10, and recorded the scaled errordecrease for the initial and final solute fraction estimates for each design (in Figure 5.30).Compared to coupling with DC resistivity data, in the heterogeneous Case 2 we obtainedon average lower errors for the final solute estimates, in the homogeneous Case 1 though,the error of initial solute estimates were higher. In both cases the error decrease for thefinal estimate was fairly similar regardless the GW design.92Figure 5.29: ADMM estimates. Upper left and right: Estimates for initial and final solutedistribution for Case 1. Bottom left and right: Estimates for initial and final solutedistribution for Case 2. Isosurfaces at ω = 0.25, 0.5 and 0.75 are displayed. The truemodels are plotted in Figure 5.35.3.2 SummaryADMM method applied to joint inversion with EM data lead to similar conclusions as whenapplied to seawater intrusion example in section 5.1 with DC resistivity data. The jointinversion with ADMM decreased the estimate errors in all cases compared to a simpler cou-pled approach only or a separate inversion. In the heterogeneous case, the ADMM approachmanaged to improve the estimates consistently even for different designs of groundwaterdata locations.The main difference compared to coupling with DC resistivity data is thus increasing93Case 1 Case 2Figure 5.30: The relative error decrease for different GW sampling designs; the dotted lineis for the ω0 relative error, the full line for the ωf relative error decrease. The plottedresults are based on different transects of wells plotted in Figure 5.10.the computational time due to the TDEM forward and inverse model. On the other hand,during our tests and implementation of the ADMM for EM data we could observe that thejoint inversion is converging to a minimum, though possibly in more iterations, regardlessthe initial guess or particular weights setup. Coupling with TDEM data therefore seemedto be more robust, compared to coupling with DC resistivity data.94Chapter 6ConclusionsJoint hydrogeophysical inversion is an attractive option to improve parameter and stateestimates for groundwater models and to better monitor and predict solute transport pro-cesses. However, when evaluating the hydrological and geophysical data, simpler coupledapproaches which use existing codes for groundwater or geophysical models are often pre-ferred. In this thesis we developed both models in the same computational environmentand derived the sensitivities analytically to reduce the computational burden of solving theinverse problem. With such a setup we were able to test different joint inversion schemesfor problems in 3D, which was one of the goals of this thesis.In chapter two, we introduced the groundwater models developed for this study. Thefirst one enables solute transport calculations, such as salt tracer tests, and the second oneis used to model more complex processes such as seawater intrusions by enabling variabledensity flow. For geophysical modeling we implemented DC resistivity and time domainelectromagnetic surveys, which are both sensitive to changes in subsurface electrical con-ductivity (and thus also fluid conductivity). DC resistivity is a popular technique amonghydrogeologists, since, compared to some other geophysical methods, it represents a rela-tively cheap method with easy data collection.In most coupled approaches, the groundwater and geophysical properties are linked viaan empirical petrophysical relationship, whose parameters need to be determined. For lab-oratory experiments, e.g. sand box with homogeneous soil, the parameter values can beeasily calibrated, but in the field conditions the parameters are uncertain at best. There-fore, when coupling the groundwater and geophysical data, we need to consider two differentsituations: first, when the petrophysical relationship is known and second, when its param-95eters are unknown and we only assume some structure similarity between solute contentand electrical conductivity. In both cases though, we minimize an objective function whichcontains the two data misfit terms, regularization terms and the coupling term comingfrom the petrophysical relationship or a structural constraint. The minimization can becomputationally challenging due to the different speeds of convergence for the groundwa-ter and geophysical models. In this work we propose computationally efficient strategiesto minimize this objective function, which we tested for multiple synthetic examples, andcompared with more standard approaches.For the joint inversion with a known petrophysical constraint, we proposed here analternating direction method of multipliers (ADMM), which splits the minimization intotwo subproblems; the GW and geophysical part. This is advantageous as the two differentdata misfits do not need to be weighted in one objective function, and we can still proceedwith the joint minimization. We tested the ADMM approach on a seawater intrusion ex-ample, where we estimated the initial and final saltwater distribution, while assuming theknowledge of other groundwater model parameters. In addition to ADMM, we examined asimpler coupled approach, where the result of a groundwater data inversion served as thereference model for the geophysical inverse problem. ADMM reduced the estimate errorsby roughly a half for the final solute estimates and to about 60% for the initial solutecompared to a simpler coupled approach errors. Similar quantitative results were obtainedwhen coupling the GW model with both DC resistivity and TDEM data. Moreover, alter-ing the GW parameters used in the inversion from the true ones, still lead to significantimprovement of ADMM estimates.For the joint inversion without a known petrophysical constraint, we implemented thecross gradient field product and joint total variation as structure similarity measures andtested the joint inversion on synthetic examples with solute tracer. To minimize the ob-jective function we applied two methods; the block coordinate descent method (BCDM),which also splits the minimization onto groundwater and geophysical part, and the jointGauss-Newton (JGN) method that minimizes both models simultaneously. Both methodswere compared to a single data inversion for groundwater data.BCDM in all cases improved or at least provided equal estimates. In the case of jointGauss-Newton minimization we need to set three different weights in the objective func-tion, which is rather difficult, hence, the success in this task varies and so do the errorsof estimates. Overall, the joint Gauss-Newton method was a less robust and reliable min-96imization technique. In terms of structure similarity measures, the coupled inversion witha cross gradient field product decreased the errors of state estimates compared to a singledata inversion. However, using Joint total variation as a structure similarity measure didnot improve the error of estimates regardless the chosen minimization method. The poorperformance of Joint total variation can probably be attributed to the scarce groundwatersampling data we had in our study.The structure coupled inversion was also applied to solve the inverse problem for pa-rameter estimation, specifically to estimate the hydraulic conductivity K. We adopted thereduced parametrization via the Karhunen-Loeve expansion for the covariance matrix ofthe conductivity field, which was assumed to be lognormal. To our knowledge, no structuralinversion has been done using a parametric approach. The improvement in terms of errorwas relatively small compared to the groundwater data inversion alone; however, the jointinversion led to lower heterogeneity contrast estimates of the conductivity K, which werecloser to the true field.In summary, based on all the synthetic studies conducted, we can confirm that thejoint inversion approach leads to better estimates of states or parameters for GW models.The goodness of estimates was largely determined by the quality and amount of dataavailable, which plays a key role in the success of solving any inverse problem. Differentdata sampling locations result in differing final estimate errors. Increasing the amount ofgroundwater data led to small improvements of the joint approach estimates compared tothe coupled approach or single data inversion alone. Results also show that, if the amountof groundwater data samples was reduced, the joint approaches were still able to givereasonable estimates with lower errors, indicating that the joint approach is less sensitiveto the amount of GW data samples.However, the fact that the joint inversion improves the estimates of GW parametersand states is not surprising or new, many coupled approaches have been documented toimprove the estimates by previous research studies. The main contribution of this worktherefore lies in introducing the minimization techniques that can handle the joint inversionfor 3D inverse problems.To our knowledge, neither ADMM or BCDM have been applied to solve hydrogeophys-ical inversion problems before. In this work we implemented both methods, including theimplicit computation of sensitivities and complex time stepping techniques, and showed thatwe can solve the inverse problems in 3D in a relatively short time. Moreover, the proposed97methods can also be applied to large scale problems, since they enable parallelization.Structure coupled inversion in hydrogeophysics is still unexplored area, therefore, wepicked a rather simple scenario to test the feasibility of such an approach. The resultsare promising, with some improvement even for the estimation of heterogeneous hydraulicconductivity. For future research in this area, it would be worthwhile to test the structurecoupled inversion on more complex examples and with different types of geophysical data.However, without the knowledge of the petrophysical relationship, we have less informationabout the system and can not expect as much improvement as in the situation where avalid petrophysical relationship is implemented.98Bibliography[1] E. Abarca, J. Carrera, X. Sa´nchez-Vila, and M. Dentz. Anisotropic dispersive henryproblem. Advances in Water Resources, 30(4):913–926, 2007. → pages x, 26[2] P. Ackerer, A. Youne`s, and M. Mancip. A new coupling algorithm for density-drivenflow in porous media. Geophysical research letters, 31(12), 2004. → pages 23[3] G. E. Archie et al. The electrical resistivity log as an aid in determining somereservoir characteristics. Transactions of the AIME, 146(01):54–62, 1942. → pages 42[4] P. M. Barlow and E. G. Reichard. Saltwater intrusion in coastal regions of northamerica. Hydrogeology Journal, 18(1):247–260, 2010. → pages 57[5] P. Bauer-Gottwein, B. N. Gondwe, L. Christiansen, D. Herckenrath, L. Kgotlhang,and S. Zimmermann. Hydrogeophysical exploration of three-dimensional salinityanomalies with the time-domain electromagnetic method (tdem). Journal ofHydrology, 380(3):318–329, 2010. → pages 6[6] J. Bear and A. H.-D. Cheng. Modeling groundwater flow and contaminant transport,volume 23. Springer Science & Business Media, 2010. → pages 16, 21[7] J. Bear, A. H.-D. Cheng, S. Sorek, D. Ouazar, and I. Herrera. Seawater intrusion incoastal aquifers: concepts, methods and practices, volume 14. Springer Science &Business Media, 1999. → pages 15[8] J. Beaujean, F. Nguyen, A. Kemna, A. Antonsson, and P. Engesgaard. Calibration ofseawater intrusion models: Inverse parameter estimation using surface electricalresistivity tomography and borehole data. Water Resources Research, 50(8):6828–6849, 2014. → pages 6, 8, 58[9] A. Beck and L. Tetruashvili. On the convergence of block coordinate descent typemethods. SIAM Journal on Optimization, 23(4):2037–2060, 2013. → pages 54[10] P. Blomgren and T. F. Chan. Color tv: total variation methods for restoration ofvector-valued images. IEEE transactions on image processing, 7(3):304–309, 1998. →pages 4599[11] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimizationand statistical learning via the alternating direction method of multipliers.Foundations and Trends R© in Machine Learning, 3(1):1–122, 2011. → pages 12, 52[12] M. B. Cardenas and M. R. Kanarek. Soil moisture variation and dynamics across awildfire burn boundary in a loblolly pine (pinus taeda) forest. Journal of Hydrology,519:490–502, 2014. → pages 3[13] M. Cardiff and P. Kitanidis. Bayesian inversion for facies detection: An extensiblelevel set framework. Water Resources Research, 45(10), 2009. → pages 10, 44, 73[14] A. Cardona, J. Carrillo-Rivera, R. Huizar-Alvarez, and E. Graniel-Castro.Salinization in coastal aquifers of arid zones: an example from santo domingo, bajacalifornia sur, mexico. Environmental Geology, 45(3):350–366, 2004. → pages 4[15] J. Carrera, J. J. Hidalgo, L. J. Slooten, and E. Va´zquez-Sun˜e´. Computational andconceptual issues in the calibration of seawater intrusion models. Hydrogeologyjournal, 18(1):131–145, 2010. → pages 2[16] M. A. Celia, T. F. Russell, I. Herrera, and R. E. Ewing. An eulerian-lagrangianlocalized adjoint method for the advection-diffusion equation. Advances in WaterResources, 13(4):187–206, 1990. → pages 19[17] M. Commer and G. A. Newman. Three-dimensional controlled-sourceelectromagnetic and magnetotelluric joint inversion. Geophysical JournalInternational, 178(3):1305–1316, 2009. → pages 50[18] M. Commer, M. B. Kowalsky, J. Doetsch, G. A. Newman, and S. Finsterle.Mpitough2: A parallel parameter estimation framework for hydrological andhydrogeophysical applications. Computers & Geosciences, 65:127–135, 2014. →pages 9[19] J.-C. Comte and O. Banton. Cross-validation of geo-electrical and hydrogeologicalmodels to evaluate seawater intrusion in coastal aquifers. Geophysical researchletters, 34(10), 2007. → pages 57[20] J.-C. Comte, O. Banton, J.-L. Join, and G. Cabioch. Evaluation of effectivegroundwater recharge of freshwater lens in small islands by the combined modeling ofgeoelectrical data and water heads. Water Resources Research, 46(6), 2010. → pages3[21] S. Costabel and U. Yaramanci. Estimation of water retention parameters fromnuclear magnetic resonance relaxation time distributions. Water resources research,49(4):2068–2079, 2013. → pages 4100[22] D. L. de Castro and R. M. G. C. Branco. 4-d ground penetrating radar monitoring ofa hydrocarbon leakage site in fortaleza (brazil) during its remediation process: a casehistory. Journal of Applied Geophysics, 54(1):127–144, 2003. → pages 4[23] W. Deng and W. Yin. On the global and linear convergence of the generalizedalternating direction method of multipliers. Journal of Scientific Computing, 66(3):889–916, 2016. → pages 52[24] J. Doetsch, N. Linde, M. Pessognelli, A. G. Green, and T. Gu¨nther. Constraining 3-delectrical resistance tomography with gpr reflection data for improved aquifercharacterization. Journal of Applied Geophysics, 78:68–76, 2012. → pages 4[25] J. E. Doherty, R. J. Hunt, and M. J. Tonkin. Approaches to highly parameterizedinversion: A guide to using PEST for model-parameter and predictive-uncertaintyanalysis. US Department of the Interior, US Geological Survey, 2011. → pages 7[26] J. Douglas, Jr and T. F. Russell. Numerical methods for convection-dominateddiffusion problems based on combining the method of characteristics with finiteelement or finite difference procedures. SIAM Journal on Numerical Analysis, 19(5):871–885, 1982. → pages 19[27] M. Droske and M. Rumpf. A variational approach to nonrigid morphological imageregistration. SIAM Journal on Applied Mathematics, 64(2):668–687, 2004. → pages44[28] M. J. Ehrhardt, K. Thielemans, L. Pizarro, D. Atkinson, S. Ourselin, B. F. Hutton,and S. R. Arridge. Joint reconstruction of pet-mri by exploiting structural similarity.Inverse Problems, 31(1):015001, 2014. → pages 45[29] G. O. Essink. Modeling three-dimensional density dependent groundwater flow at theisland of texel, the netherlands. Coastal Aquifer Management-Monitoring, Modeling,and Case Studies, page 77, 2016. → pages 57[30] C. G. Farquharson and D. W. Oldenburg. A comparison of automatic techniques forestimating the regularization parameter in non-linear inverse problems. GeophysicalJournal International, 156(3):411–425, 2004. → pages 48[31] D. V. Fitterman. Mapping saltwater intrusion in the biscayne aquifer, miami-dadecounty, florida using transient electromagnetic sounding. Journal of Environmentaland Engineering Geophysics, 19(1):33–43, 2014. → pages 57[32] D. E. Fowler and S. M. Moysey. Estimation of aquifer transport parameters fromresistivity monitoring data within a coupled inversion framework. Journal ofHydrology, 409(1):545–554, 2011. → pages 3101[33] L. A. Gallardo and M. A. Meju. Characterization of heterogeneous near-surfacematerials by joint 2d inversion of dc resistivity and seismic data. GeophysicalResearch Letters, 30(13), 2003. → pages 45[34] L. A. Gallardo and M. A. Meju. Joint two-dimensional dc resistivity and seismictravel time inversion with cross-gradients constraints. Journal of GeophysicalResearch: Solid Earth, 109(B3), 2004. → pages 10[35] L. A. Gallardo and M. A. Meju. Structure-coupled multiphysics imaging ingeophysical sciences. Reviews of Geophysics, 49(1), 2011. → pages 10, 44[36] E. Ghadimi, A. Teixeira, I. Shames, and M. Johansson. Optimal parameter selectionfor the alternating direction method of multipliers (admm): quadratic problems.IEEE Transactions on Automatic Control, 60(3):644–658, 2015. → pages 12, 52[37] E. Haber. Computational methods in geophysical electromagnetics, volume 1. SIAM,2014. → pages 18, 30, 33, 48[38] E. Haber and M. H. Gazit. Model fusion and joint inversion. Surveys in Geophysics,34(5):675–695, 2013. → pages 44, 45[39] E. Haber and D. Oldenburg. Joint inversion: a structural approach. Inverseproblems, 13(1):63, 1997. → pages 10, 44, 45[40] F. Haeni. Application of seismic-refraction techniques to hydrologic studies. USGovernment Printing Office, 1988. → pages 4[41] P. C. Hansen. Discrete inverse problems: insight and algorithms, volume 7. Siam,2010. → pages 62[42] J. D. Hem. Study and interpretation of the chemical characteristics of natural water,volume 2254. Department of the Interior, US Geological Survey, 1985. → pages 43[43] H. R. Henry. Interfaces between salt water and fresh water in coastal aquifers. USGeological Survey Water-Supply Paper, pages C35–70, 1964. → pages 24[44] D. Herckenrath, G. Fiandaca, E. Auken, and P. Bauer-Gottwein. Sequential and jointhydrogeophysical inversion using a field-scale groundwater model with ert and tdemdata. Hydrology and Earth System Sciences, 17(10):4043–4060, 2013. → pages 3, 6, 7[45] D. Herckenrath, N. Odlum, V. Nenna, R. Knight, E. Auken, and P. Bauer-Gottwein.Calibrating a salt water intrusion model with time-domain electromagnetic data.Groundwater, 51(3):385–397, 2013. → pages 3, 6, 9, 58102[46] T. Hermans, A. Vandenbohede, L. Lebbe, R. Martin, A. Kemna, J. Beaujean, andF. Nguyen. Imaging artificial salt water infiltration using electrical resistivitytomography constrained by geostatistical data. Journal of Hydrology, 438:168–180,2012. → pages 7[47] A. Hinnell, T. Ferre´, J. Vrugt, J. Huisman, S. Moysey, J. Rings, and M. Kowalsky.Improved extraction of hydrologic information from geophysical data through coupledhydrogeophysical inversion. Water resources research, 46(4), 2010. → pages 3, 6, 7[48] W. Hu, A. Abubakar, and T. M. Habashy. Joint electromagnetic and seismicinversion using structural constraints. Geophysics, 74(6):R99–R109, 2009. → pages 10[49] J. Irving and K. Singha. Stochastic inversion of tracer test and electrical geophysicaldata to estimate hydraulic conductivities. Water Resources Research, 46(11), 2010.→ pages 3, 6, 8[50] A. Jardani, A. Revil, A. Bole`ve, A. Crespy, J.-P. Dupont, W. Barrash, andB. Malama. Tomography of the darcy velocity from self-potential measurements.Geophysical Research Letters, 34(24), 2007. → pages 42[51] A. Jardani, A. Revil, and J. Dupont. Stochastic joint inversion of hydrogeophysicaldata for salt tracer test monitoring and hydraulic conductivity imaging. Advances inWater Resources, 52:62–77, 2013. → pages 3, 7, 8, 80[52] D. H. Jayawickreme, C. S. Santoni, J. H. Kim, E. G. Jobba´gy, and R. B. Jackson.Changes in hydrology and salinity accompanying a century of agricultural conversionin argentina. Ecological Applications, 21(7):2367–2379, 2011. → pages 3[53] D. H. Jayawickreme, E. G. Jobba´gy, and R. B. Jackson. Geophysical subsurfaceimaging for ecological applications. New Phytologist, 201(4):1170–1175, 2014. →pages 3[54] T. C. Johnson, R. J. Versteeg, H. Huang, and P. S. Routh. Data-domain correlationapproach for joint hydrogeologic inversion of time-lapse hydrogeologic andgeophysical data. Geophysics, 74(6):F127–F140, 2009. → pages 10, 44, 73[55] A. Kacimov, M. Sherif, J. Perret, and A. Al-Mushikhi. Control of sea-water intrusionby salt-water pumping: Coast of oman. Hydrogeology Journal, 17(3):541–558, 2009.→ pages 57[56] A. Kemna, A. Binley, F. Day-Lewis, A. Englert, B. Tezkan, J. Vanderborght,H. Vereecken, and P. Winship. Solute transport processes. In AppliedHydrogeophysics, pages 117–159. Springer, 2006. → pages 42, 43103[57] O. Kolditz, R. Ratke, H.-J. G. Diersch, and W. Zielke. Coupled groundwater flowand transport: 1. verification of variable density flow and transport models.Advances in Water Resources, 21(1):27–46, 1998. → pages 15, 22[58] L. F. Konikow, D. J. Goode, and G. Z. Hornberger. A three-dimensionalmethod-of-characteristics solute-transport model (MOC3D). US Geological Survey,1996. → pages 19[59] M. B. Kowalsky, S. Finsterle, J. Peterson, S. Hubbard, Y. Rubin, E. Majer, A. Ward,and G. Gee. Estimation of field-scale soil hydraulic and dielectric parameters throughjoint inversion of gpr and hydrological data. Water Resources Research, 41(11), 2005.→ pages 41[60] C. D. Langevin and W. Guo. Modflow/mt3dms–based simulation of variable-densityground water flow and transport. Ground Water, 44(3):339–351, 2006. → pages 19[61] C. D. Langevin and M. Zygnerski. Effect of sea-level rise on salt water intrusion neara coastal well field in southeastern florida. Groundwater, 51(5):781–803, 2013. →pages 57[62] C. D. Langevin, D. T. Thorne Jr, A. M. Dausman, M. C. Sukop, and W. Guo.Seawat version 4: A computer program for simulation of multi-species solute andheat transport. Technical report, Geological Survey (US), 2008. → pages 19[63] P. G. Lelie`vre and D. W. Oldenburg. A comprehensive study of including structuralorientation information in geophysical inversions. Geophysical Journal International,178(2):623–637, 2009. → pages 10[64] P. G. Lelie`vre, D. W. Oldenburg, and N. C. Williams. Integrating geological andgeophysical data through advanced constrained inversions*. Exploration Geophysics,40(4):334–341, 2009. → pages 10[65] T. Lochbu¨hler, J. Doetsch, R. Brauchler, and N. Linde. Structure-coupled jointinversion of geophysical and hydrological data. Geophysics, 78(3):ID1–ID14, 2013. →pages 11, 45, 73[66] J. Mead. Parameter estimation: A new approach to weighting a priori information.J. Inv. Ill-posed Problems, 16(2):175–194, 2008. → pages 33[67] A. Medina and J. Carrera. Coupled estimation of flow and solute transportparameters. Water Resources Research, 32(10):3063–3076, 1996. → pages 48[68] T. Mills, P. Hoekstra, M. Blohm, and L. Evans. Time domain electromagneticsoundings for mapping sea-water intrusion in monterey county, california. GroundWater, 26(6):771–782, 1988. → pages 57104[69] D. M. Molodtsov, V. N. Troyan, Y. V. Roslov, and A. Zerilli. Joint inversion ofseismic traveltimes and magnetotelluric data with a directed structural constraint.Geophysical Prospecting, 61(6):1218–1228, 2013. → pages 10[70] M. Moorkamp. Integrated Imaging of the Earth: Theory and Applications, volume218. John Wiley & Sons, 2016. → pages 44[71] L. Neilson-Welch and L. Smith. Saline water intrusion adjacent to the fraser river,richmond, british columbia. Canadian geotechnical journal, 38(1):67–82, 2001. →pages xi, 59, 60[72] V. Nenna, D. Herckenrath, R. Knight, N. Odlum, and D. McPhee. Application andevaluation of electromagnetic methods for imaging saltwater intrusion in coastalaquifers: Seaside groundwater basin, california. Geophysics, 78(2):B77–B88, 2013. →pages 4, 57[73] W. Nijland, M. Van der Meijde, E. A. Addink, and S. M. De Jong. Detection of soilmoisture and vegetation water abstraction in a mediterranean natural area usingelectrical resistivity tomography. Catena, 81(3):209–216, 2010. → pages x, 3, 4[74] R. Nishihara, L. Lessard, B. Recht, A. Packard, and M. I. Jordan. A general analysisof the convergence of admm. arXiv preprint, 2015. → pages 53[75] J. Nocedal and S. Wright. Numerical optimization. Springer Science & BusinessMedia, 2006. → pages 36, 51, 54[76] D. S. Oliver and Y. Chen. Recent progress on reservoir history matching: a review.Computational Geosciences, 15(1):185–221, 2011. → pages 80[77] J. G. Paine. Determining salinization extent, identifying salinity sources, andestimating chloride mass using surface, borehole, and airborne electromagneticinduction methods. Water Resources Research, 39(3), 2003. → pages 4[78] E. E. Poeter, M. C. Hill, E. R. Banta, S. Mehl, and S. Christensen. Ucode 2005 andsix other computer codes for universal sensitivity analysis, calibration, anduncertainty evaluation constructed using the jupiter api. Technical report, 2006. →pages 7[79] D. Pollock and O. A. Cirpka. Fully coupled hydrogeophysical inversion of syntheticsalt tracer experiments. Water Resources Research, 46(7), 2010. → pages 8, 80[80] D. Pollock and O. A. Cirpka. Fully coupled hydrogeophysical inversion of alaboratory salt tracer experiment monitored by electrical resistivity tomography.Water Resources Research, 48(1), 2012. → pages 7, 8105[81] M. Putti and C. Paniconi. Picard and newton linearization for the coupled model forsaltwater intrusion in aquifers. Advances in water resources, 18(3):159–170, 1995. →pages 24[82] B. S. RamaRao, A. M. LaVenue, G. De Marsily, and M. G. Marietta. Pilot pointmethodology for automated calibration of an ensemble of conditionally simulatedtransmissivity fields: 1. theory and computational experiments. Water ResourcesResearch, 31(3):475–493, 1995. → pages 80[83] M. Razaviyayn, M. Hong, and Z.-Q. Luo. A unified convergence analysis of blocksuccessive minimization methods for nonsmooth optimization. SIAM Journal onOptimization, 23(2):1126–1153, 2013. → pages 54[84] Y. Rubin and S. S. Hubbard. Hydrogeophysics, volume 50. Springer Science &Business Media, 2006. → pages 3[85] T. F. Russell and M. A. Celia. An overview of research on eulerian–lagrangianlocalized adjoint methods (ellam). Advances in Water Resources, 25(8):1215–1231,2002. → pages 19[86] E. Sanchez-Leo´n, C. Leven, C. Haslauer, and O. Cirpka. Combining 3d hydraulictomography with tracer tests for improved transport characterization. Groundwater,2015. → pages 80[87] W. E. Sanford and J. P. Pope. Current challenges using models to forecast seawaterintrusion: lessons from the eastern shore of virginia, usa. Hydrogeology Journal, 18(1):73–93, 2010. → pages 57[88] G. Sapiro and D. L. Ringach. Anisotropic diffusion of multivalued images withapplications to color filtering. IEEE transactions on image processing, 5(11):1582–1586, 1996. → pages 45[89] S. Sorek and V. Borisov. Modified eulerian–lagrangian formulation for hydrodynamicmodeling. Journal of Computational Physics, 231(8):3083–3100, 2012. → pages 19[90] K. Steklova and E. Haber. 3d joint hydrogeophysical inversion using similaritymeasures. submitted to Applied Mathematics and Computations. → pages iv[91] K. Steklova and E. Haber. Joint hydrogeophysical inversion: state estimation forseawater intrusion models in 3d. Computational Geosciences, pages 1–20, 2016. →pages iv, 13[92] D. K. Todd and L. W. Mays. Groundwater hydrology edition, 1980. → pages 16106[93] F. Trabelsi, A. B. Mammou, J. Tarhouni, C. Piga, and G. Ranieri. Delineation ofsaltwater intrusion zones using the time domain electromagnetic method: thenabeul–hammamet coastal aquifer case study (ne tunisia). Hydrological Processes, 27(14):2004–2020, 2013. → pages 4, 57[94] C. R. Vogel. Computational methods for inverse problems, volume 23. Siam, 2002. →pages 62[95] C. I. Voss, C. T. Simmons, and N. I. Robinson. Three-dimensional benchmark forvariable-density flow and transport simulation: matching semi-analytic stabilitymodes for steady unstable convection in an inclined porous box. Hydrogeologyjournal, 18(1):5–23, 2010. → pages 21[96] J.-M. Vouillamoz, S. Sokheng, O. Bruyere, D. Caron, and L. Arnout. Towards abetter estimate of storage properties of aquifer with magnetic resonance sounding.Journal of Hydrology, 458:51–58, 2012. → pages 4[97] K. Vozoff and D. Jupp. Joint inversion of geophysical data. Geophysical JournalInternational, 42(3):977–991, 1975. → pages 41[98] A. S. Ward, M. N. Gooseff, M. Fitzgerald, T. J. Voltz, and K. Singha. Spatiallydistributed characterization of hyporheic solute transport during baseflow recessionin a headwater mountain stream using electrical geophysical imaging. Journal ofHydrology, 517:362–377, 2014. → pages 3[99] B. Wohlberg, D. M. Tartakovsky, and M. Dentz. Linearized functional minimizationfor inverse modeling. In XIX International Conference on Water Resources(CMWR), Urbana-Champaign, IL, USA, 2012. → pages 12107
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computational methods in hydrogeophysics
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computational methods in hydrogeophysics Steklova, Klara 2017
pdf
Page Metadata
Item Metadata
Title | Computational methods in hydrogeophysics |
Creator |
Steklova, Klara |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | Parameter and state estimation for groundwater models within a coupled hydrogeophysical framework has become common in the last few years as it has been shown that such estimates are usually better than those from a single data inversion. Different approaches have been suggested in literature to combine the essentially two different modalities in order to obtain better estimates for groundwater models, and improve monitoring of processes such as solute transport. However, the coupled approaches usually come at a price of higher computational cost and difficulties in coupling the geophysical and groundwater inverse problems. Unlike in other studies, we developed both the groundwater and geophysical models in the same computational environment in order to test different minimization strategies. When solving the coupled inverse problem, the objective function consists of data misfit and regularization terms as well as a coupling term that relates groundwater and geophysical states. We present a novel approach to solve the inverse problem using an Alternating Direction Method of Multipliers (ADMM) to minimize the coupled objective function. ADMM enables us to treat the groundwater and geophysical part separately and thus use existing software with minor changes. However, ADMM as well as many other coupled approaches relies on implementing some petrophysical relationship to couple the groundwater and geophysical variable. Such relationships are usually uncertain and hard to parametrize for a large region and can potentially produce solute mass errors in the final model estimates. Therefore, in this thesis we examine coupled approaches that replace the fixed petrophysical relationship by a more loose structure similarity constraint. Besides, we propose efficient computational methods to minimize the objective function when there is no explicit petrophysical constraint. All approaches were tested on 3D synthetic examples. In the solute tracer test we estimated hydraulic conductivity or solute distribution using a structure coupled inversion, and were able to reduce the errors compared to a single data inversion alone. For a more complex example of seawater intrusion we implemented the ADMM method, and obtained better estimates for the solute distribution compared to just considering each data separately, or solving the problem with a simple coupled approach. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-03-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0343090 |
URI | http://hdl.handle.net/2429/60815 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2017-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2017_may_steklova_klara.pdf [ 4.07MB ]
- Metadata
- JSON: 24-1.0343090.json
- JSON-LD: 24-1.0343090-ld.json
- RDF/XML (Pretty): 24-1.0343090-rdf.xml
- RDF/JSON: 24-1.0343090-rdf.json
- Turtle: 24-1.0343090-turtle.txt
- N-Triples: 24-1.0343090-rdf-ntriples.txt
- Original Record: 24-1.0343090-source.json
- Full Text
- 24-1.0343090-fulltext.txt
- Citation
- 24-1.0343090.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0343090/manifest