strange error for 20 km event depth

[lpk-py]

I'm extracting seismograms from a customed axisem database using instaseis. It works good, but each time I set event depth at 20 km, I got this error

     47         s = db.get_seismograms(src, rec, components='Z', kind='displacement', \
---> 48                   remove_source_shift=False, return_obspy_stream=False, dt=dt)

/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_instaseis_db.pyc in get_seismograms(self, source, receiver, components, kind, remove_source_shift, reconvolve_stf, return_obspy_stream, dt, kernelwidth)
    246         # Call the _get_seismograms() method of the respective implementation.
    247         data = self._get_seismograms(source=source, receiver=receiver,
--> 248                                      components=components)
    249
    250         if dt is None:

/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_netcdf_instaseis_db.pyc in _get_seismograms(self, source, receiver, components)
    201         coordinates = Coordinates(s=rotmesh_s, phi=rotmesh_phi, z=rotmesh_z)
    202
--> 203         element_info = self._get_element_info(coordinates=coordinates)
    204
    205         return self._get_data(

/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_netcdf_instaseis_db.pyc in _get_element_info(self, coordinates)
    126                     break
    127             else:  # pragma: no cover
--> 128                 raise ValueError("Element not found")
    129
    130             if not self.read_on_demand:

ValueError: Element not found

[martinvandriel]

This error seems to suggest you request outside of the domain of where you stored the wavefields, but without more detailed information (input files on the axisem run) it's hard to say how this happens.

[lpk-py]

I saved the whole kernel wavefields in axisem, from surface to center of the core.

I extracted seismograms for sources at varying depth. Any depth (e.g. 0, 10, 30, 50, 400, 1200, 3400, 6000 km) works, except 20 km (up to now no others were found problematic).

# minimal and maximal colatitude in degree for kernel wavefields
# (only for dumptype displ_only)
KERNEL_COLAT_MIN   0.
KERNEL_COLAT_MAX   180.

# minimal and maximal radius in km for kernel wavefields
# (only for dumptype displ_only)
KERNEL_RMIN           0.  
KERNEL_RMAX        7000.

inparam_advanced.txt

[martinvandriel]

This is a bit strange, is it exactly 20 only, i.e. would 20.1 or 19.9 work? And which model do you use (inparam_mesh, inparam_basic would help)?

[martinvandriel]

Is it the same for all distances as well?

[lpk-py]

Sorry, I broke my conda env after updating some packages.

I'll post the results after.

[lpk-py]

I tested depths from 16 to 25 km. 16, 17, 18, 23, 24, 25 km are ok. So problems are with from 18.x to 22.x km.

It is both related to source depth and distance. I tried 0:30:180 deg for evdp=19km. The error raises at dis=90deg.

get seismograms of 7 receivers from source@19.0km
dis=0.0
dis=30.0
dis=60.0
dis=90.0
---------------------------------------------------------------------------
......
/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_instaseis_db.pyc in get_seismograms(self, source, receiver, components, kind, remove_source_shift, reconvolve_stf, return_obspy_stream, dt, kernelwidth)
    246         # Call the _get_seismograms() method of the respective implementation.
    247         data = self._get_seismograms(source=source, receiver=receiver,
--> 248                                      components=components)
    249
    250         if dt is None:

/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_netcdf_instaseis_db.pyc in _get_seismograms(self, source, receiver, components)
    201         coordinates = Coordinates(s=rotmesh_s, phi=rotmesh_phi, z=rotmesh_z)
    202
--> 203         element_info = self._get_element_info(coordinates=coordinates)
    204
    205         return self._get_data(

/home/user/anaconda2/lib/python2.7/site-packages/instaseis/database_interfaces/base_netcdf_instaseis_db.pyc in _get_element_info(self, coordinates)
    126                     break
    127             else:  # pragma: no cover
--> 128                 raise ValueError("Element not found")
    129
    130             if not self.read_on_demand:

ValueError: Element not found

[lpk-py]

The mesh model is ak135f, with a period of 16s and 20 cores.

I didn't keep the mesh files. But you can generate easily with these parameters.

[lpk-py]

I tested evdp=18:0.1:19km, dis=0:10:180deg and evdp=23:-0.1:22deg, dis=0:10:180deg.

The error raises from evdp=18.2km, dis=20deg and evdp=22.4km, dis=90deg, respectively.

[lpk-py]

18.2 to 22.4 km depth, so the error raises when the source is just beneath the Moho in ak135f. Is it related to meshes at these depths?

[martinvandriel]

I suspect that in this this layer with doubling the correct element is not necessarily in the first six closest midpoints. Increasing the value to 8 seems to fix it.

https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_netcdf_instaseis_db.py#L74

[martinvandriel]

you will need to do the developers installation to use the patch until the next release

[krischer]

I'd increase this even a bit more. In the given case I can see how 8 might be needed - but even a sligtht distortion might require even more, especially for example with very thin crustal layer. Unlikely on Earth but maybe happens on other planets? I cannot see any harm is just increasing this to something much larger, for example 64.