[CIG-SHORT] RELAX Question

[Austin Madson]

Hi Sylvain,


Thanks for the response!

The lake in question has seen a marked increase in water level in the last
several years - I expect to see between ~2-3cm of downward defo on the
ground from DInSAR (and other geodetic studies). However, most studies
using models have seen a modeled elastic response 2-6x lower than the above
value (This is for the example load I have been using [+2.44m of water].

The elastic model I was using recently (REAR
<https://github.com/danielemelini/rear>), which computes Green's functions
from load deformation coefficients along the lines of Farrel (1972), gave
outputs along the lines of ~6x too low of an elastic response from the same
load. However, that model does not account for viscoelastic responses and I
was hoping to use RELAX to model that.

That's why I was expecting much lower values from the elastic response. I
figured they would be different, but not nearly that much different!

In response to keeping things in meters and MPa:
If my Lamé params are 30000 (MPa), gamma is 8.33e-7 (meters) and my load is
2.39e-2 (MPa) [which is equivalent to ~2.44m of water load] - the response
is attached below. Is this value in meters as is, or do we divide the
surface defo by 1000 to get meters?

Either way, the response is at least 10x too great (if the attached output
is in meters) or at least 20x too little (if we divide the attached output
by 1000 and it's in meters). Based on what I know from other models and
studies, I expect 2-6x too low of the modeled elastic response from this
load.

Maybe I'm missing something somewhere (a decimal place is off, or
something)?

I hope this makes sense.

Cheers,
Austin M

On Fri, Jan 29, 2016 at 10:50 PM, Sylvain Barbot <sylbar.vainbot at gmail.com>
wrote:

> Hi Austin,
>
> You can call me Sylvain. I think I am going to bounce it back to you and
> ask you how you arrived at different expectations for the displacement. The
> displacement is a function of the load and the area where it is applied.
> And to be sure, try to run your model using SI units. It should remove some
> confusion. Use meters for all distance units and MPa for all stress and
> load units.
>
> Cheers,
> Sylvain
>
>
> On Saturday, January 30, 2016, Austin Madson <amadson at ucla.edu> wrote:
>
>> Prof. Barbot,
>>
>> Ahh, yes. Good catch! For some reason I was under the impression that x1
>> and x2 were the upper left coords and the length and width were the length
>> in width away from the upper left coords. Thanks for the catch!
>>
>> However, after fixing that overlap, the response was only reduced by 3cm
>> (and still well over what I was expecting to see here). Any ideas as to why
>> it's so high?
>>
>> Cheers,
>> Austin M
>>
>> On Wed, Jan 27, 2016 at 11:55 PM, Sylvain Barbot <
>> sylbar.vainbot at gmail.com> wrote:
>>
>>> Hi Austin,
>>>
>>> Your input file is set correctly. (But recall that length is in the
>>> north direction and width is in the east direction.) Relax computes
>>> the elastic response if the viscosity is infinite. From intuition
>>> alone, I am not shocked by your results. But perhaps you did not
>>> intend to overlap your patches? Right now they overlap for 5km.
>>>
>>> Cheers,
>>> Sylvain
>>>
>>> On Thu, Jan 28, 2016 at 2:43 PM, Austin Madson <amadson at ucla.edu> wrote:
>>> > Prof. Barbot,
>>> >
>>> > Thanks for the well written explanation!
>>> >
>>> > Before I move on to viscous layers, I ran a sanity check (see below)
>>> to get
>>> > the response from a 2.44m water level increase in a lake. The response
>>> was
>>> > much higher than anticipated (~19cm where it should be <2cm). Is RELAX
>>> not
>>> > set up to properly look at elastic outputs (Which is fine, because I'm
>>> after
>>> > viscoelastic outputs). But, if the initial elastic response is this
>>> high I'm
>>> > concerned the viscoelastic response will still be too great (with a
>>> tau = ~
>>> > 10-13 years)). I'm quite certain everything is correct here. What do
>>> you
>>> > think? I have attached a PDF from the output as well.
>>> >
>>> > Thanks in advance!
>>> >
>>> > Also, do you know anyone using RELAX for surface loading?
>>> >
>>> > Cheers,
>>> > Austin M
>>> >
>>> > Lamé params in GPa, load in MPa, sampling in km.
>>> >
>>> > # grid dimension (sx1,sx2,sx3)
>>> > 512 512 512
>>> > # sampling (dx1,dx2,dx3), smoothing (beta, nyquist)
>>> > 0.9765625 0.9765625 0.9765625 0.2 0
>>> > # origin position (x0,y0) and rotation
>>> > 0 0 0
>>> > # observation depth (displacement and stress)
>>> > 0 0
>>> > # output directory
>>> > $WDIR
>>> > # lambda, mu, gamma (gamma = (1 - nu) rho g / mu)
>>> > 30 30 8.33e-4
>>> > # time interval, (positive time step) or (negative skip, scaling)
>>> > 0 1 1
>>> >
>>> > # number of observation planes
>>> > 0
>>> > # number of observation points
>>> > 0
>>> > # number of stress observation segments
>>> > 0
>>> > # number of prestress interfaces
>>> > 0
>>> > # number of linear viscous interfaces
>>> > 0
>>> > # number of nonlinear viscous interfaces
>>> > 0
>>> > # number of fault creep interfaces
>>> > 0
>>> > # number of inter-seismic strike-slip segments
>>> > 0
>>> > # number of inter-seismic tensile segments
>>> > 0
>>> > # number of events
>>> > 1
>>> > # number of coseismic strike-slip segments
>>> > 0
>>> > # number of coseismic tensile segments
>>> > 0
>>> > # number of coseismic dilatation point sources
>>> > 0
>>> > # number of surface loads
>>> > #890730
>>> > 2
>>> > # nb x1 x2 length(x) width(y) t3 T phi (t3 in MPa)
>>> > 1 15 -15 30 30 0.023912  0 0
>>> > 2 10 -35 20 20 0.023912  0 0
>>> > EOF
>>> >
>>> >
>>> >
>>> >
>>> > On Tue, Jan 26, 2016 at 3:07 PM, Sylvain Barbot <
>>> sylbar.vainbot at gmail.com>
>>> > wrote:
>>> >>
>>> >> Hi Austin,
>>> >>
>>> >> If your spatial sampling is in km then your Lamé parameters are in
>>> GPa,
>>> >> your load in MPa and your displacement in meters. If your spatial
>>> sampling
>>> >> was in meters, then your Lamé parameters and your load would be in
>>> GPa (or
>>> >> more precisely, the same, arbitrary, unit). The difference comes from
>>> stress
>>> >> being rigidity times strain, and strain being displacement divided by
>>> >> distance. If the distance and displacements are not in the same
>>> units, you
>>> >> bias the strain.
>>> >>
>>> >> Cheers,
>>> >> Sylvain
>>> >>
>>> >>
>>> >> On Wednesday, January 27, 2016, Austin Madson <amadson at ucla.edu>
>>> wrote:
>>> >>>
>>> >>> Prof. Barbot,
>>> >>>
>>> >>> Thanks again for the response. Just to confirm, if in the below
>>> example
>>> >>> (the same as before) the Lamé params are in GPa (30 and 30) then t3
>>> would be
>>> >>> in GPa (of the load) and the output surface displacement would still
>>> be in
>>> >>> meters?
>>> >>>
>>> >>> Thanks!!
>>> >>>
>>> >>> Cheers,
>>> >>> Austin M
>>> >>>
>>> >>>
>>> >>> # grid dimension (sx1,sx2,sx3)
>>> >>> 512 512 512
>>> >>> # sampling (dx1,dx2,dx3), smoothing (beta, nyquist)
>>> >>> 0.9765625 0.9765625 0.9765625 0.2 0
>>> >>> # origin position (x0,y0) and rotation
>>> >>> 0 0 0
>>> >>> # observation depth (displacement and stress)
>>> >>> 0 0
>>> >>> # output directory
>>> >>> $WDIR
>>> >>> # lambda, mu, gamma (gamma = (1 - nu) rho g / mu)
>>> >>> 30 30 8.33e-4
>>> >>> # time interval, (positive time step) or (negative skip, scaling)
>>> >>> 1 1 1
>>> >>> # number of observation planes
>>> >>> 0
>>> >>> # number of observation points
>>> >>> 0
>>> >>> # number of stress observation segments
>>> >>> 0
>>> >>> # number of prestress interfaces
>>> >>> 0
>>> >>> # number of linear viscous interfaces
>>> >>> 0
>>> >>> # number of nonlinear viscous interfaces
>>> >>> 0
>>> >>> # number of fault creep interfaces
>>> >>> 0
>>> >>> # number of inter-seismic strike-slip segments
>>> >>> 0
>>> >>> # number of inter-seismic tensile segments
>>> >>> 0
>>> >>> # number of events
>>> >>> 1
>>> >>> # number of coseismic strike-slip segments
>>> >>> 0
>>> >>> # number of coseismic tensile segments
>>> >>> 0
>>> >>> # number of coseismic dilatation point sources
>>> >>> 0
>>> >>> # number of surface loads
>>> >>> 1
>>> >>> # nb x1 x2 length width t3 T phi
>>> >>> 1 0 0 3.937 3.937 20 0 0
>>> >>> EOF
>>> >>>
>>> >>> On Fri, Jan 22, 2016 at 4:30 AM, Sylvain Barbot
>>> >>> <sylbar.vainbot at gmail.com> wrote:
>>> >>>>
>>> >>>> Hi Austin,
>>> >>>>
>>> >>>> With your input file, you compute the displacement field at the
>>> >>>> surface with units of meters, due to a 4x4 km patch where a load of
>>> 20
>>> >>>> MPa is applied. Your spatial sampling is 4 km. The value of gamma
>>> >>>> should be 8.33e-4 in units of 1/km.
>>> >>>>
>>> >>>> Cheers,
>>> >>>> Sylvain
>>> >>>>
>>> >>>> On Fri, Jan 22, 2016 at 4:24 PM, Austin Madson <amadson at ucla.edu>
>>> wrote:
>>> >>>> > Prof. Barbot,
>>> >>>> >
>>> >>>> > Thanks again for the response - you've been  very helpful thusfar.
>>> >>>> >
>>> >>>> > Just to confirm:
>>> >>>> > Trying to get at displacement units (let's call them km). If my
>>> inputs
>>> >>>> > are
>>> >>>> > the following, and I divide out the 1e3 from the displacement
>>> field -
>>> >>>> > the
>>> >>>> > units are just that, km. Is that correct? Or, if I keep the Lamé
>>> >>>> > parameters
>>> >>>> > (along with t3)  in GPa instead of kPa - does the Gamma value
>>> then =
>>> >>>> > 8.33e2
>>> >>>> > [gamma=8.33e2 = (1-0.25)*3400*9.8/30 or gamma=8.33e4 =
>>> >>>> > (1-0.25)*3400*9.8/30000000 ],  and we still divide out the 1e3
>>> from
>>> >>>> > the
>>> >>>> > displacement to get the units in "km"?
>>> >>>> >
>>> >>>> >
>>> >>>> >
>>> >>>> > Cheers and thanks!
>>> >>>> >
>>> >>>> > Austin M
>>> >>>> >
>>> >>>> >
>>> >>>> >
>>> >>>> >
>>> >>>> > # grid dimension (sx1,sx2,sx3)
>>> >>>> > 512 512 512
>>> >>>> > # sampling (dx1,dx2,dx3), smoothing (beta, nyquist)
>>> >>>> > 0.9765625 0.9765625 0.9765625 0.2 0
>>> >>>> > # origin position (x0,y0) and rotation
>>> >>>> > 0 0 0
>>> >>>> > # observation depth (displacement and stress)
>>> >>>> > 0 0
>>> >>>> > # output directory
>>> >>>> > $WDIR
>>> >>>> > # lambda, mu, gamma (gamma = (1 - nu) rho g / mu)
>>> >>>> > 30 30 8.33e-4
>>> >>>> > # time interval, (positive time step) or (negative skip, scaling)
>>> >>>> > 1 1 1
>>> >>>> > # number of observation planes
>>> >>>> > 0
>>> >>>> > # number of observation points
>>> >>>> > 0
>>> >>>> > # number of stress observation segments
>>> >>>> > 0
>>> >>>> > # number of prestress interfaces
>>> >>>> > 0
>>> >>>> > # number of linear viscous interfaces
>>> >>>> > 0
>>> >>>> > # number of nonlinear viscous interfaces
>>> >>>> > 0
>>> >>>> > # number of fault creep interfaces
>>> >>>> > 0
>>> >>>> > # number of inter-seismic strike-slip segments
>>> >>>> > 0
>>> >>>> > # number of inter-seismic tensile segments
>>> >>>> > 0
>>> >>>> > # number of events
>>> >>>> > 1
>>> >>>> > # number of coseismic strike-slip segments
>>> >>>> > 0
>>> >>>> > # number of coseismic tensile segments
>>> >>>> > 0
>>> >>>> > # number of coseismic dilatation point sources
>>> >>>> > 0
>>> >>>> > # number of surface loads
>>> >>>> > 1
>>> >>>> > # nb x1 x2 length width t3 T phi
>>> >>>> > 1 0 0 3.937 3.937 20 0 0
>>> >>>> > EOF
>>> >>>> >
>>> >>>> >
>>> >>>> > On Tue, Jan 19, 2016 at 11:20 PM, Sylvain Barbot
>>> >>>> > <sylbar.vainbot at gmail.com>
>>> >>>> > wrote:
>>> >>>> >>
>>> >>>> >> Hi Austin,
>>> >>>> >>
>>> >>>> >> I often get questions about units. There are no formal choice of
>>> >>>> >> units
>>> >>>> >> in Relax, except that S.I. units, or more precisely, their
>>> >>>> >> inter-relationships, is assumed. The expected input of the
>>> current
>>> >>>> >> version of Relax is
>>> >>>> >>
>>> >>>> >> n x1 x2 length width t3 T phi
>>> >>>> >>
>>> >>>> >> where n is the index, x1 and x2 and north and east, length and
>>> width
>>> >>>> >> are the dimension of the source patch, t3 is the traction (force
>>> per
>>> >>>> >> unit area), and T and phi is the period and phase assuming the
>>> >>>> >> following model
>>> >>>> >>
>>> >>>> >> p(t)=t3*sin(2*pi*t/period+phi)
>>> >>>> >>
>>> >>>> >> t3 is the load in units of stress. If you prescribed your Lamé
>>> >>>> >> parameters in units of MPa, you need to describe the load in MPa.
>>> >>>> >>
>>> >>>> >> I apologise that the examples I sent you correspond to an old
>>> version
>>> >>>> >> of Relax. The current version assumes stress units. (The previous
>>> >>>> >> version expected t3/G.) The factor of 1,000 comes in depending on
>>> >>>> >> your
>>> >>>> >> choice of length units. If you use km, it introduces a factor of
>>> >>>> >> 1,000
>>> >>>> >> in the displacement field. You can compensate this effect by
>>> dividing
>>> >>>> >> the rigidity by the same factor.
>>> >>>> >>
>>> >>>> >> Cheers,
>>> >>>> >> Sylvain
>>> >>>> >>
>>> >>>> >>
>>> >>>> >> On Wed, Jan 20, 2016 at 2:03 PM, Austin Madson <amadson at ucla.edu
>>> >
>>> >>>> >> wrote:
>>> >>>> >> > Prof. Barbot,
>>> >>>> >> >
>>> >>>> >> > Thanks for the response and the example data. I am working out
>>> my
>>> >>>> >> > parameters
>>> >>>> >> > (and units now). I think I have a lock on the elastic and
>>> >>>> >> > viscoelastic
>>> >>>> >> > params. However, I am having trouble with the load parameters.
>>> >>>> >> >
>>> >>>> >> > I see that on some of your input files the surface
>>> load/traction
>>> >>>> >> > parameters
>>> >>>> >> > are;
>>> >>>> >> > nb x1 x2 length width t3 T phi
>>> >>>> >> >
>>> >>>> >> > and another has;
>>> >>>> >> > nb x1 x2 t3 length width T phi
>>> >>>> >> >
>>> >>>> >> > I assume the former is correct, yes?
>>> >>>> >> >
>>> >>>> >> > Also, I'm trying to wrap my head around the units for
>>> force(t3). In
>>> >>>> >> > the
>>> >>>> >> > GRACE example (the _km.xyz file), the z is used in the
>>> following
>>> >>>> >> > calculation
>>> >>>> >> > ($3*9.8/3e10) in order to get force(t3). What are the initial
>>> units
>>> >>>> >> > in
>>> >>>> >> > the
>>> >>>> >> > _km.xyz file for the third column, z? And what are the units
>>> after
>>> >>>> >> > the
>>> >>>> >> > calculation ($3*9.8/3e10) in order to derive force(t3)?
>>> >>>> >> >
>>> >>>> >> > Also, it appears that the output vertical (z) displacement
>>> units
>>> >>>> >> > are in
>>> >>>> >> > mm,
>>> >>>> >> > may I assume that is correct?
>>> >>>> >> >
>>> >>>> >> > Cheers and thanks A LOT for your responses!
>>> >>>> >> > Austin M
>>> >>>> >> >
>>> >>>> >> > On Tue, Jan 19, 2016 at 9:29 AM, Sylvain Barbot
>>> >>>> >> > <sylbar.vainbot at gmail.com>
>>> >>>> >> > wrote:
>>> >>>> >> >>
>>> >>>> >> >> Hi Austin,
>>> >>>> >> >>
>>> >>>> >> >> The Tibet lakes have a great loading signal. You can use
>>> Relax to
>>> >>>> >> >> simulate the deformation, but you need to be careful about how
>>> >>>> >> >> important the elastic stratification is. If your data are
>>> close to
>>> >>>> >> >> the
>>> >>>> >> >> source, you're fine. I presume you are using paleo-shoreline
>>> data.
>>> >>>> >> >> If
>>> >>>> >> >> so, use the shallow rigidity for the entire domain. If you are
>>> >>>> >> >> modeling far-field data, you should use a code that have a
>>> layered
>>> >>>> >> >> elastic structure, i.e., (Farrell, 1972).
>>> >>>> >> >>
>>> >>>> >> >> The example attached models the deformation from the Indian
>>> >>>> >> >> monsoon
>>> >>>> >> >> using Grace data. The example cycleload.sh creates a periodic
>>> >>>> >> >> surface
>>> >>>> >> >> load.
>>> >>>> >> >>
>>> >>>> >> >> Best wishes,
>>> >>>> >> >> Sylvain
>>> >>>> >> >>
>>> >>>> >> >>
>>> >>>> >> >>
>>> >>>> >> >> On Tue, Jan 19, 2016 at 10:51 AM, Austin Madson <
>>> amadson at ucla.edu>
>>> >>>> >> >> wrote:
>>> >>>> >> >> > Professor Barbot,
>>> >>>> >> >> >
>>> >>>> >> >> > Thanks for the prompt response. I looked at the
>>> documentation
>>> >>>> >> >> > and the
>>> >>>> >> >> > man
>>> >>>> >> >> > page (as well as the examples on the PDF on the geodynamics
>>> >>>> >> >> > website)
>>> >>>> >> >> > and
>>> >>>> >> >> > have made progress.
>>> >>>> >> >> >
>>> >>>> >> >> > As an fyi, we plan on modeling several hundred years of
>>> surface
>>> >>>> >> >> > loading
>>> >>>> >> >> > from
>>> >>>> >> >> > very large lakes (i.e. time dependent loads). Are there any
>>> >>>> >> >> > RELAX
>>> >>>> >> >> > surface
>>> >>>> >> >> > load examples laying around (they can just be
>>> single/multiple
>>> >>>> >> >> > point
>>> >>>> >> >> > and
>>> >>>> >> >> > even
>>> >>>> >> >> > instantaneous (no time function)? The GRACE example on the
>>> >>>> >> >> > example
>>> >>>> >> >> > pdf
>>> >>>> >> >> > "hides" the data in various dat files (which makes sense if
>>> >>>> >> >> > you're
>>> >>>> >> >> > using
>>> >>>> >> >> > a
>>> >>>> >> >> > lot of data, but it's hard to see their example inputs).
>>> >>>> >> >> >
>>> >>>> >> >> > I expect to see only millimeters of deformation throughout
>>> the
>>> >>>> >> >> > run -
>>> >>>> >> >> > will
>>> >>>> >> >> > the outputs be precise enough if I use the correct SI units
>>> >>>> >> >> > throughout?
>>> >>>> >> >> >
>>> >>>> >> >> > Cheers,
>>> >>>> >> >> > Austin Madson
>>> >>>> >> >> >
>>> >>>> >> >> > On Sun, Jan 17, 2016 at 4:06 PM, Sylvain Barbot
>>> >>>> >> >> > <sylbar.vainbot at gmail.com>
>>> >>>> >> >> > wrote:
>>> >>>> >> >> >>
>>> >>>> >> >> >> Hi Austin,
>>> >>>> >> >> >>
>>> >>>> >> >> >> The details of the input file for surface loads are in the
>>> >>>> >> >> >> documentation
>>> >>>> >> >> >> and in the man page. There are some pdf files on the
>>> >>>> >> >> >> geodynamics
>>> >>>> >> >> >> website
>>> >>>> >> >> >> that show examples of surface loads for the Himalayan
>>> region if
>>> >>>> >> >> >> I
>>> >>>> >> >> >> recall
>>> >>>> >> >> >> well.
>>> >>>> >> >> >>
>>> >>>> >> >> >> But the input file simply expects a list of squares with
>>> their
>>> >>>> >> >> >> associated
>>> >>>> >> >> >> traction. So depending on your project, this can be a
>>> single
>>> >>>> >> >> >> point,
>>> >>>> >> >> >> say
>>> >>>> >> >> >> to
>>> >>>> >> >> >> represent the loading of a dam, or a complicated function
>>> of
>>> >>>> >> >> >> space
>>> >>>> >> >> >> and
>>> >>>> >> >> >> time.
>>> >>>> >> >> >> If your change of load is instantaneous, you need only one
>>> >>>> >> >> >> event. Of
>>> >>>> >> >> >> your
>>> >>>> >> >> >> load is time dependent, you need as many events as you have
>>> >>>> >> >> >> time
>>> >>>> >> >> >> steps.
>>> >>>> >> >> >> I
>>> >>>> >> >> >> can be more specific, but I need more information about
>>> your
>>> >>>> >> >> >> practical
>>> >>>> >> >> >> goals.
>>> >>>> >> >> >>
>>> >>>> >> >> >> Cheers,
>>> >>>> >> >> >> Sylvain
>>> >>>> >> >> >>
>>> >>>> >> >> >>
>>> >>>> >> >> >> On Monday, January 18, 2016, Austin Madson <
>>> amadson at ucla.edu>
>>> >>>> >> >> >> wrote:
>>> >>>> >> >> >>>
>>> >>>> >> >> >>> Professor Barbot,
>>> >>>> >> >> >>>
>>> >>>> >> >> >>> I have a quick question re: RELAX - Is there any more
>>> >>>> >> >> >>> information
>>> >>>> >> >> >>> out
>>> >>>> >> >> >>> there with respect to utilizing RELAX for deformation
>>> >>>> >> >> >>> responses to
>>> >>>> >> >> >>> surface
>>> >>>> >> >> >>> loading? I have looked through all of the examples as
>>> well as
>>> >>>> >> >> >>> the
>>> >>>> >> >> >>> slides/tutorials on the geodynamics.org website.
>>> >>>> >> >> >>>
>>> >>>> >> >> >>> Can you provide any further information? Or can you point
>>> me
>>> >>>> >> >> >>> in a
>>> >>>> >> >> >>> better
>>> >>>> >> >> >>> direction?
>>> >>>> >> >> >>>
>>> >>>> >> >> >>> Cheers,
>>> >>>> >> >> >>> Austin Madson
>>> >>>> >> >> >
>>> >>>> >> >> >
>>> >>>> >> >
>>> >>>> >> >
>>> >>>> >
>>> >>>> >
>>> >>>
>>> >>>
>>> >
>>>
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20160130/790f9ed8/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 000-up.grd-plot.pdf
Type: application/pdf
Size: 40779 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20160130/790f9ed8/attachment-0001.pdf>