Rayleigh Waves Power Spectrum of Vertical Displacement¶

In [1]:

Copied!





## Distributed python packages
import os
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

from math import radians, log
from sys import exit

from wmsan.subfunctions_rayleigh_waves import loop_SDF, download_ww3_local, open_bathy

__author__ = "Lisa Tomasetto"
__copyright__ = "Copyright 2024, UGA"
__credits__ = ["Lisa Tomasetto"]
__version__ = "0.1"
__maintainer__ = "Lisa Tomasetto"
__email__ = "lisa.tomasetto@univ-grenoble-alpes.fr"
__status__ = "Production"
## Distributed python packages
import os
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

from math import radians, log
from sys import exit

from wmsan.subfunctions_rayleigh_waves import loop_SDF, download_ww3_local, open_bathy

__author__ = "Lisa Tomasetto"
__copyright__ = "Copyright 2024, UGA"
__credits__ = ["Lisa Tomasetto"]
__version__ = "0.1"
__maintainer__ = "Lisa Tomasetto"
__email__ = "lisa.tomasetto@univ-grenoble-alpes.fr"
__status__ = "Production"

Rayleigh Source of microseisms¶

This program aims at modelizing the ambient noise source in the secondary microseismic range for Rayleigh waves. It is based on Ardhuin et al. (2011) and Stutzmann et al. (2012) articles. We will compute the equivalent source for the power spectrum of the vertical displacement $S_{DF} \text{ in m.s}$.

$$S_{DF}(f_s) \approx \frac{2\pi f_s C}{\rho_s^2 \beta^5} F_{p3D}(k_2 \approx 0, f_s)$$

where:

$f_s$ is the seismic frequency in Hz (twice the ocean wave frequency)
C is the amplitude response functions for the normal modes, from Longuet-Higgins (1950). $C = \sum_{i=1}^{4} c_i^2$, no dimension.
$\rho_s$ is the rock density of the crust. Here $\rho_s = 2830 \text{ kg.m}^{-3}$
$\beta$ is the shear wave velocity. Here $\beta = 2.8 \text{ km.s}^{-1}$
$F_{p3D}(k_2 \approx 0, f_s)$ the spectral density of the pressure field at the ocean surface or directional wave spectra in $ \text{Pa}^2.\text{m}^2.\text{s}$.

In our case Rayleigh waves are dominant in the periodic band from 5s to 12s so we will integrate this over the corresponding frequency band.

$$S_{DF} = \int_{0.083}^{0.200} S_{DF}(f_s) df_s$$

The values of the Rayleigh waves site effect coefficient can be found in Longuet-Higgins (1950), the values are tabulated in longuet_higgins.txt file. The $F_{p3D}$ can be retrieved using oceanic waves model resources from the Ifremer institute in Brest, France. We provide a cell downloading the files directly.

References:

Parameters¶

This implementation is depending on a few parameters. Some are described above in the formula of Rayleigh wave noise source but others are for the user to choose.

Physical Constants¶

Default values are given below. Please change with your own values if needed.

In [2]:

Copied!





# seismic waves
VS_CRUST= 2.8  # km/s
RHO_S = 2830  # kg/m3
F1 = 1/12  # frequency to integrate from
F2 = 1/2  # frequency to integrate to
parameters = [VS_CRUST, RHO_S, F1, F2]
# seismic waves
VS_CRUST= 2.8  # km/s
RHO_S = 2830  # kg/m3
F1 = 1/12  # frequency to integrate from
F2 = 1/2  # frequency to integrate to
parameters = [VS_CRUST, RHO_S, F1, F2]

Dates¶

Then the dates the user wants to focus on, loops on hours, days, months and years are available setting an empty bracket symbol '[ ]'.

In [3]:

Copied!





# dates
YEAR = 2008  # year of interest
MONTH = [3]  # loop if array, compute all months of the year if empty list []
DAY = np.arange(8, 20)  # loop if array, compute all days of the month if empty list []
HOUR = []  # loop if array, compute every 3 hours of the day if empty list []

date_vec = [YEAR, MONTH, DAY, HOUR]
# dates
YEAR = 2008  # year of interest
MONTH = [3]  # loop if array, compute all months of the year if empty list []
DAY = np.arange(8, 20)  # loop if array, compute all days of the month if empty list []
HOUR = []  # loop if array, compute every 3 hours of the day if empty list []

date_vec = [YEAR, MONTH, DAY, HOUR]

Spatial Extent¶

In [4]:

Copied!





# extent
lat_min = -78 # -78 min
lat_max = 80 # 80 max
lon_min = -180 # -180 min
lon_max = 180 # 180 max

extent = [lon_min, lon_max, lat_min, lat_max]
# extent
lat_min = -78 # -78 min
lat_max = 80 # 80 max
lon_min = -180 # -180 min
lon_max = 180 # 180 max

extent = [lon_min, lon_max, lat_min, lat_max]

Additional Files¶

See detailed information in cells below.

Paths to Files¶

In [9]:

Copied!





# ftp path of WW3 data
ftp_path_to_files = "ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/SISMO/GLOBAL05_%d_REF102040/"%YEAR

# local path for WW3 data
ww3_local_path = "../../data/%d"%YEAR  # path where the data will be downloaded

# Longuet-Higgins site effect coefficients
longuet_higgins_file = '../../data/longuet_higgins.txt'

# bathymetry default
file_bathy = "../../data/LOPS_WW3-GLOB-30M_dataref_dpt.nc"  #0.5x0.5 degree grid bathymetry

paths = [file_bathy, ww3_local_path, longuet_higgins_file]
# ftp path of WW3 data
ftp_path_to_files = "ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/SISMO/GLOBAL05_%d_REF102040/"%YEAR

# local path for WW3 data
ww3_local_path = "../../data/%d"%YEAR  # path where the data will be downloaded

# Longuet-Higgins site effect coefficients
longuet_higgins_file = '../../data/longuet_higgins.txt'

# bathymetry default
file_bathy = "../../data/LOPS_WW3-GLOB-30M_dataref_dpt.nc"  #0.5x0.5 degree grid bathymetry

paths = [file_bathy, ww3_local_path, longuet_higgins_file]

Download WW3 Files¶

For the model files go to ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/SISMO/, then choose the year(s) and month(s) corresponding files. We download the directional wave spectra file, extension p2l.nc (default). It will download a .nc file containing the full output of WW3 (including significant waveheight Hs).

To download other versions add the prefix option, for example: prefix = 'SWOT_WW3-GLOB-30M'.

In [10]:

Copied!

download_ww3_local(YEAR, MONTH, ftp_path_to_files, ww3_local_path)
download_ww3_local(YEAR, MONTH, ftp_path_to_files, ww3_local_path)

[Errno 17] File exists: '../../data/2008'
Downloading can take some time...

-----------------------------------------------------------------

ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/SISMO/GLOBAL05_2008_REF102040/WW3-GLOB-30M_200803_p2l.nc downloaded

-----------------------------------------------------------------

WW3 files downloaded in ../../data/2008
current directory :  /Users/tomasetl/GITLAB/ww3-source-maps/notebooks/rayleigh_waves

2024-07-02 11:23:09 URL: ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/SISMO/GLOBAL05_2008_REF102040/WW3-GLOB-30M_200803_p2l.nc [1404618101] -> "WW3-GLOB-30M_200803_p2l.nc" [1]

Bathymetry file¶

The bathymetry (or waterlevel) is necessary to compute the site effect for a given phase.

Two bathymetry grids are available for this notebook:

(default) "../../data/LOPS_WW3-GLOB-30M_dataref_dpt.nc": a 0.5°x0.5° bathymetry file corresponding to WW3 hindcast resolution.
(to download) a 1 arcmin resolution ETOPOv2 bathymetry netcdf file.(refined_bathymetry = True)

Both file should be located in the ww3-source-maps/data/ directory.

ETOPOv2 file is also available here: https://www.ngdc.noaa.gov/thredds/catalog/global/ETOPO2022/60s/60s_bed_elev_netcdf/catalog.html?dataset=globalDatasetScan/ETOPO2022/60s/60s_bed_elev_netcdf/ETOPO_2022_v1_60s_N90W180_bed.nc

[WARNING] use this refined bathymetry on small grids otherwise memory errors might occur (typically 30° lat x 30° lon)

If you wish to use your own bathymetry file:

latitude in ° 1D array should be named zlat.
longitude in ° 1D array should be named zlon.
depth in meters 2D grid should be named dpt1.

In [11]:

Copied!

dpt1, zlon, zlat = open_bathy(file_bathy, refined_bathymetry=False, extent=extent)
dpt1, zlon, zlat = open_bathy(file_bathy, refined_bathymetry=False, extent=extent)

Save and Plot¶

If you want to save the 3-hourly matrix as a netcdf file set the save variable to True. save = True

If you want to plot the maps and save them as .png files set the plot_type variable to:

plot_type = 'hourly' for plots every 3-hours (WW3 resolution)
plot_type = 'daily' for summed daily plots
plot_type = 'monthly' for summed monthly plots
plot_type = 'yearly' for summed yearly plots

In [12]:

Copied!

save = True  # save matrix as netcdf
plot_type = 'hourly' # plot Force 'hourly', 'daily', 'monthly', 'yearly'
save = True  # save matrix as netcdf
plot_type = 'hourly' # plot Force 'hourly', 'daily', 'monthly', 'yearly'

Main Loop¶

The next cell computes ambient noise Rayleigh waves source maps depending on the previously defined parameters.

To use other WW3 files than default add the prefix option, for example: prefix = 'SWOT_WW3-GLOB-30M'.

In [13]:

Copied!





loop_SDF(paths, dpt1, zlon, zlat,
         date_vec = date_vec,
         extent = extent,
         parameters=parameters,
         plot_type=plot_type,
         save = save, vmin=0, vmax=1e-17)
loop_SDF(paths, dpt1, zlon, zlat,
         date_vec = date_vec,
         extent = extent,
         parameters=parameters,
         plot_type=plot_type,
         save = save, vmin=0, vmax=1e-17)

File WW3  ../../data/2008/WW3-GLOB-30M_200803_p2l.nc
Rayleigh source maps done!