University of Alaska Fairbanks
Department of Geosciences

GEOS F493 / F693 - Geodetic Methods

Lab 13: Volcano Deformation

Note that you DO NOT have to work on LuNGS today!

Introduction

Today you will use some co-eruptive displacements from Redoubt volcano's 2009 eruption and hopefully learn a bit about its source. You will do 2 things: (1) implement the Mogi model, (2) find the best source parameters given the displacements.

Caveats:

(1) I've already shown that the Mogi source isn't the best model for the co-eruptive phase. We'll use it anyway for the data, because it's a simple model, and we're here to learn the methods.

(2) The method I'll have you implement to solve for the non-linear parts is inefficient and imprecise. Again, it's one way to gain an intuition to find solutions to these problems; you should embark on a journey to learn about the many efficient non-linear methods / solvers that are out there yourself; often they take a parameter range in which they do their magic, which reduces to a single function call. Using this wouldn't result in much learning in a lab.

I bring up these points to make sure you don't assume the methods used here are the ones you should use for the rest of your career; far from it! They're a starting point.

The Forward Model

First, start out with implementing a function [dr, dz] = mogi(x,y,d,C), which - given a set of distances (x,y) in carthesian coordinates away from a source location, a source depth d and a source strength C - returns radial and horizontal displacements at these points according to the Mogi model (check slides for formulas).

You can test your model with these values:

>> [dr, dz] = mogi([1000, -1000, 2000, 0]', [1000, -1000, 2000, 0]', 17000, (0.75*0.05*1000^3)/pi)

dr =

   0.003400610640462
   0.003400610640462
   0.006596168094920
                   0


dz =

   0.040878112348779
   0.040878112348779
   0.039645659112996
   0.041303185923502			
		

The x,y distances are in meters away from a source at 17 km depth that changes in volume by 0.05 km3. What are the units of dr, dz?

The ``Inverse'' Model

Download the file invert_mogi_incomplete.m. It contains a matrix displacements; its contents is described in the comments and a vector volcano which contains in UTM coordinates (meters) the location of the volcano.

The network setup is shown in the figure here:

Write a script that searches within +/-10 km of the volcano's location, within a depth range of, say, 40 km for a Mogi source. Note that the Mogi model is non-linear with respect to location and hence you'll need to come up with a non-linear solution scheme. I'd recommend you do a grid search in 1 km intervals through your model domain. You can invert for your source strength C using least-squares since the model is linear in that parameter.

Note that this is not the best way to solve this problem, but it provides considerable intuition as to what more sophisticated methods might be after.

Hints

A few things you need to do (Test each step individually! Do not try to slay the entire beast in one fell swoop!):

• setup your data vector concatenating ur, uz in correct units
• turn the global coordinates of the stations to radial distances with respect to the volcano location (subtract volcano location)
• set a new variable current_resid_norm=inf
• write loops to iterate over x, y, depth
• set up your design matrix (G) to solve for C in the innermost loop (depends on r, d)
• invert for C - the source strength
• compute the forward model using the mogi function you've written before. The parameter values for this function call will be the x,y,d values from the triple loop you are in. C will come from the least squares solution just calculated.
• calculate the norm of the residual vector (use function norm in Matlab)
• if your new norm is better than any of the previous norms, save your current model parameters in variables, say, best_x, best_y, best_d, best_C; save this currently best norm in current_resid_norm

What is the expected source depth and the volume change of the magma reservoir (convert from source strength to volume (see lecture slides and above)!)? Give these values in sensible units.

Deliverables: (submit via blackboard in a zip file!)

• mogi function file
• inverse_mogi script
• separate document with answers to questions in bold (PDF!)