Model times kernel is not travel time. about mc_kernel HOT 30 CLOSED

A basic question regarding the units of the kernels we compute: Our kernels are essentially given by eqn. 5.13 in Tarjes thesis, right? So their unit should be s^2 / m . When multiplied with the background model what should come out is the absolute traveltime. But is this consistent with eqn. 5.12 in Tarjes thesis which relates differential traveltimes with (percentual) model perturbations via (relative?) kernels. I probably don't quite understand the difference between absolute and relative kernels, can sb explain?

Here are some more basic things that can/should probably be checked to make sure that our kernels are plausible (not all directly related to this issue reported by Simon):

• Are the kernel values (after dividing by the element volumes) independent of which mesh is used?
• In which range do the kernel values lie? Are they similar to published kernels? In the literature the kernels are usually plotted on scales b/w -5 to 5 * 10-6 s/km3
• Are the kernel values independent of the source magnitude used in the fwd simulation?
• Are there unrealistically large values close to the source and the receiver, which may obscure the integrated kernel values? Looking at the K_x array it seems that there are some elements with extremely large values, which may completely dominate the summation with the background model.
• Do the kernels look the same for more "complicated" (i.e. ones where the source is not at the northpole) source-receiver geometries?
• Do the kernels visually look similar to published kernels?
• Do the kernel computations work for all possible phases?

from mc_kernel.

Just some quick thoughts:

Are the kernel values (after dividing by the element volumes) independent of which mesh is used?

• I think so! at least by assuming a very fine initial mesh, it should be possible to project it onto other inversion grids. I am saying very fine because otherwise we have not sampled the domain correctly.

In which range do the kernel values lie? Are they similar to published kernels? In the literature the kernels are usually plotted on scales b/w -5 to 5 * 10-6 s/km3

• I assume that the best is to check the integral of that (as we have already talked about) and compare it with the travel time...but I have also seen the same values. They actually report it like that to be independent of the inversion grid (if I am not mistaken).

Are the kernel values independent of the source magnitude used in the fwd simulation?

• Yes! This is why it is divided by the length of Vp (last term in eq 5.13)

from mc_kernel.

@auerl:

I think in 5.12 and 5.13 the perturbations are absolute, as this comes from 5.3. To understand the integral with the model, it is maybe useful to think about the analogy with rays. However, I have not seen a rigorous proof and I could not make it up myself in 5 minutes. When you come across it, please let me know.

@kasra-hosseini:

Kernel values: I guess the difference is, whether we talk about the kernel as a function of space K(x) or the values Aij in the matrix. The first ones should be independent of the mesh, the latter depend on the normalization chosen such that it can be multiplied with the modelvector.

from mc_kernel.

@martinvandriel
You are right, I was just thinking in this way:

\int_{Volume}KdV ---> that should be in second (arrival time) which is A (as you said) in dT = Ac

If we only calculate K, it should be then sec/km3; However, we dont do this. So if we go for fine inversion grid and then divide it by volume of each (sort of similar to numerical integration), we should be able to say it is actually K (ad hoc), and can be projected on any other inversion grid (which is relatively larger than the initial one). Does that make sense?

from mc_kernel.

I think to just look at the kernel values it would be easier to just use an orthonormal basis like Ludwigs voxels, because in this case the relation between the kernel and the matrix is trivial, see the draft of the kernel paper, eq 7. We could then just plot the kernels with the tools we have and look at the values.

from mc_kernel.

Agree! I remember that derivation from Warnemuende...that would be indeed the best.

from mc_kernel.

@auerl

Are there unrealistically large values close to the source and the receiver, which may obscure the
integrated kernel values? Looking at the K_x array it seems that there are some elements with
extremely large values, which may completely dominate the summation with the background
model.

Can you quantify that? The values at source and receiver should be considerably larger, since all energy has to go through these cells, which gives them a large influence on the result.

@martinvandriel

I think to just look at the kernel values it would be easier to just use an orthonormal basis like
Ludwigs voxels, because in this case the relation between the kernel and the matrix is trivial, see
the draft of the kernel paper, eq 7. We could then just plot the kernels with the tools we have and
look at the values.

I agree. The code can calculate volumetric kernels on any mesh, so we might also do that with tetrahedrons. But either should work.

from mc_kernel.

@auerl
I might as well answer all your questions. The answers are about the current state of the code (3c93755):

Are the kernel values (after dividing by the element volumes) independent of which mesh is used?

Difficult to say. Should test that

In which range do the kernel values lie? Are they similar to published kernels? In the literature the kernels are usually plotted on scales b/w -5 to 5 * 10-6 s/km3

See Martin's answer on Kx vs Aij. The only plots I know are of Kx.

Are the kernel values independent of the source magnitude used in the fwd simulation?

They should be and as I see it, they are.

Do the kernels look the same for more "complicated" (i.e. ones where the source is not at the northpole) source-receiver geometries?

Yes, that's fine.

Do the kernels visually look similar to published kernels?

Yes, they do

Do the kernel computations work for all possible phases?

I tested it a bit, but that could be done more rigorously. It did not really work well for more complicated phases with the time-domain integration for
$\int K_X(t) dt$. Could be tested with the Parseval-theorem based integration in frequency domain.

So to summarize @auerl, @kasra-hosseini :
Tasks for the inclined would be:

• Test more complicated phases
• Test different meshes. To remove the multiplication by volume, set the volume term in l. 367 in slave.f90 https://github.com/sstaehler/kerner/blob/master/slave.f90#L367 to 1:

    call int_kernel(ibasisfunc)%initialize_montecarlo(parameters%nkernel, &
    1, &
    parameters%allowed_error, &
    parameters%allowed_relative_error) 

from mc_kernel.

We (Ludwig, Dr. Boehm, me) tried to write down why the kernel times the model should be absolute traveltime and I think we did not quite get there. Does anyone of you have notes about it? If not, do you know someone who has them and we could ask? Without understanding this properly, there is IMHO no point in using this as a test.

from mc_kernel.

Hi all,

I just compared actual kernel values computed at cscs and our machine in Zurich, and surprisingly there is quite a difference. See the figures below. No idea what this could be. The compiler is gfortran in both cases (4.8.1 in lugano and 4.8.2 in zurich), wavefields are the same, version of the code is the same, inparam file is the same, netcdf version is 4.3.1 in lugano and 4.1.3 in zurich.

Maybe an issue related to using different (incompatible) versions of netcdf to write / read the netcdf wavefields. When running the code in Zurich also the projected background model velocities are 6 (!) orders too large. From what I see, the kernel values (if not multiplied with the volume of the grid cells) does not seem to be influenced by the chosen mesh), which is good:

Another thing that I noticed is that at places where the reference model has jumps (CMB, surface, probably also the internal discontinuities) the default MC-integration parameters that are set for projecting the model onto the inversion grid don't seem to ensure convergence, since there is lateral variations. See the following screenshot:

@ALL: I think, what would be extremely helpful would be to have a running version of the Princeton's group (Dahlen) kernel code on our machine in Zurich, to which everybody has access. Would allow us to output values at various locations for similar kernels, and help to get a feeling for how large they need to be. Also, it would be good to see how exactly the test of multiplying kernels with the background model to get absolute traveltimes is implemented.

@kasra-hosseini: Could you install such a test-version for us? I think you are the only one of us who knows the code.

from mc_kernel.

@auerl
Sure, I will try to do it this weekend. Is that OK?

@martinvandriel
I am sure I saw it somewhere in Nolet's book, but unfortunately I do not have the book now. I just wrote it down and you can find the scan version here...It is for ray-based tomography.

from mc_kernel.

@kasra-hosseini
In the first equation, you only integrate the kernel itself, not the kernel times the background model? Where is the K in the following?

from mc_kernel.

Aren't the attached couple of lines based on born scattering sufficient, at least intuitively? 

-------- Original message --------
From: Martin van Driel [email protected]
Date: 12/03/2015 13:53 (GMT+00:00)
To: sstaehler/kerner [email protected]
Subject: Re: [kerner] Model times kernel is not travel time. (#21)

@kasra-hosseini
In the first equation, you only integrate the kernel itself, not the kernel times the background model? Where is the K in the following?

—
Reply to this email directly or view it on GitHub.

from mc_kernel.

@tnissen attachment?

from mc_kernel.

Attachment didn't seem to work... Here:
http://seis.earth.ox.ac.uk/20150312_141808.jpg

On 12/03/2015 14:23, Martin van Driel wrote:

@tnissen https://github.com/tnissen attachment?

—
Reply to this email directly or view it on GitHub
#21 (comment).

Tarje

<>--<>--<>--<>--<>--<>
Dept. of Earth Sciences
Oxford University
South Parks Road
Oxford OX1 3AN; UK
tel: +44 1865 282149
fax: +44 1865 272072
web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk
<>--<>--<>--<>--<>--<>

from mc_kernel.

@martinvandriel
Well the very first equation (at top) that I wrote was really only defining the question, not carefully written. However, about your second question: I think the whole kernel in ray theory is:

-\int_{ray}{1./V}ds

Not sure actually...

from mc_kernel.

@tnissen
well now we have a - b = c - d. That is obviously true if a = c and b = d, but we need to show that this is the only solution.

from mc_kernel.

Sorry, I meant: -1/V is the kernel (ray theory)...without integral.

Just some thoughts: In all sensitivity kernels: K_rho, K_lambda and K_mu we have:

-\int_{0}^{t} something dt ---> that something will be normalized and if we consider it is constant (just approximate), we will end up with: -\int_{0}^{t}dt which is -T (travel time). Does that make sense? Or I am making myself more and more confused :D.

from mc_kernel.

I think you are now confusing the integral over time and volume. The one over volume should be the relevant one here.

from mc_kernel.

You are right! similar to -\int_{ray}{1./V}ds = -T in ray theory, the relevant one should be the spatial integration.

from mc_kernel.

My line of thought at the moment is to use travelime Kernels with respect to slowness instead of velocity or moduli. For these it is straightforward to show the property we are looking for. The missing step is now to relate these to the other kernels in first order. I'll try to write it down properly...

from mc_kernel.

I would say Born scattering validates that... ?

On 12/03/2015 14:36, Martin van Driel wrote:

@tnissen https://github.com/tnissen
well now we have a - b = c - d. That is obviously true it a = c and b
= d, but we need to show that this is the only solution.

—
Reply to this email directly or view it on GitHub
#21 (comment).

Tarje

<>--<>--<>--<>--<>--<>
Dept. of Earth Sciences
Oxford University
South Parks Road
Oxford OX1 3AN; UK
tel: +44 1865 282149
fax: +44 1865 272072
web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk
<>--<>--<>--<>--<>--<>

from mc_kernel.

The more I think about it, the less clear it becomes: what is T in a finite frequency sense? delta T is clear, as it comes from comparing two signals, but what is the absolute traveltime T?

from mc_kernel.

... another thing we are struggling with right now concerns the units of our kernels: When looking at eqns. 5.12 and 5.13 in Tarjes thesis (on which our implementation is based on), there doesn't seem to be correspondence between the left and the right side in 5.12 in terms of the physical units. Can somebody explain what I am missing here? In 5.12 \delta T is given in [s] and \delta rho is given in [kg/m^3] so the kernels as defined in 5.13 ought to have the unit [s/kg]. However, if I see it correctly, 5.13 has the unit [s^2/m] (since vp is the velocity seismogram, the K inside the time integral is in [s] and the normalization factor is in [s/m^2]. Usually the kernels in all kernel galleries are plotted in the unit [s/m^3](or [km^3]).

I have the same problem in the banana donought section of Tromps paper on "Seismic tomography, adjoint methods, time reversal and banana-donought kernels" (2004, GJI). The following equations don't seem to make sense. In eqn. 46 [s] is equated with [kgm^2/s].

Does somebody have an idea?

from mc_kernel.

I will be in Zurich next week (wed afternoon/thur morning), let's sit together on this, i can't beforehand in more detail. When writing my thesis I spent quite some time exactly on this...as I recall (almost 10 years ago!) it is related to the fact that the backward/Adjoint wavefield has weird units... definitely not those of a seismic wavefield

-------- Original message --------
From: auerl [email protected]
Date: 13/03/2015 12:26 (GMT+00:00)
To: sstaehler/kerner [email protected]
Cc: tnissen [email protected]
Subject: Re: [kerner] Model times kernel is not travel time. (#21)

... another thing we are struggling with right now concerns the units of our kernels: When looking at eqns. 5.12 and 5.13 in Tarjes thesis (on which our implementation is based on), there doesn't seem to be correspondence between the left and the right side in 5.12 in terms of the physical units. Can somebody explain what I am missing here? In 5.12 \delta T is given in [s] and \delta rho is given in [kg/m^3] so the kernels as defined in 5.13 ought to have the unit [s/kg]. However, if I see it correctly, 5.13 has the unit [s^2/m] (since vp is the velocity seismogram, the K inside the time integral is in [s] and the normalization factor is in [s/m^2]. Usually the kernels in all kernel galleries are plotted in the unit s/m^3.

I have the same problem in the banana donought section of Tromps paper on "Seismic tomography, adjoint methods, time reversal and banana-donought kernels" (2004, GJI). The following equations don't seem to make sense. In eqn. 46 [s] is equated with [kgm^2/s].

Uploading trompeq.png . . .

Does somebody have an idea?

—
Reply to this email directly or view it on GitHub.

from mc_kernel.

Hi all,

It seems that what I missed before was that the delta function in the force term has a unit. See Carl Tapes thesis at page 184: "In practice, the adjoint sources have a m^−3 or m^−2 quantity as well, due to the 3D or 2D volume for the delta function, δ(x), that is applied at the source location" or wikipedia "Eine direkte Folgerung aus der Skalierungseigenschaft ist die Dimension bzw. Maßeinheit der Delta-Distribution. Sie entspricht genau der reziproken Dimension ihres Arguments".

Accounting for the additional 1/m^3 the units in the equations in Tromps paper make sense and one gets the [s^2/kg m] = [1/N] as the unit of the adjoint wavefield, if I see it correctly.

Greets,

Ludwig

from mc_kernel.

I went through the units of our implementation and everything seems consistent. The kernels, as we have them right now, are for absolute perturbations, and not relative ones. For absolute ones we just need to add an extra multiplication with the background model (see 2nd page). Should we introduce an option in the inparam file for that? Relative ones are probably more important. Having the units of all the involved quantities, I'll now check whether they are all approximately on the right order.

from mc_kernel.

Yes, that's how I had it in the original kernel code as well - it was
indeed related to the units of the backward wavefield.

The 7 orders of magnitude w.r.t. Dahlen... could that have to do with
Mij units between cgs to SI? My recollection of the Dahlen 2000 paper is
that the kernel represents an explosion (isotropic radiation), but I
don't recall the incorporation of the moment. Is the backward wavefield
for a delta function of the correct integral equal 1? At least in an
earlier version we multiplied it with a large number to avoid having all
too small numbers in the wavefield. I bet you take care of all this,
just checking.

On 17/03/2015 15:25, auerl wrote:

I went through the units of our implementation and everything seems
consistent.That the kernels, as we have them right now, are for
absolute perturbations, and not relative ones. For absolute ones we
just need to add an extra multiplication with the background model
(see 2nd page). Having the units of all the involved quantities, I'll
next check whether they are all approximately on the right order.

image1
https://cloud.githubusercontent.com/assets/5452469/6690557/84e18a72-ccc1-11e4-852b-ff3fadc6240f.jpeg

image2
https://cloud.githubusercontent.com/assets/5452469/6690559/890e7dda-ccc1-11e4-9348-229032c72631.jpeg

—
Reply to this email directly or view it on GitHub
#21 (comment).

Tarje

<>--<>--<>--<>--<>--<>
Dept. of Earth Sciences
Oxford University
South Parks Road
Oxford OX1 3AN; UK
tel: +44 1865 282149
fax: +44 1865 272072
web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk
<>--<>--<>--<>--<>--<>

from mc_kernel.

The scalar moment should not affect the kernel value, since the convolved wavefield is divided by the squared seismogram in 5.13.
I'll run a test to see whether this is the case in our code.

On the amplitude of the backward field: Good question, I'll have a look

from mc_kernel.

I close this discussion now. The underlying assumption seems to have been wrong.

from mc_kernel.