Detailed crack wave analysis¶
Here we demonstrate the crack wave analysis described in the main supplementary document. The purpose of the analysis is to determine peak positions of overtones in the spectra of crack waves recorded at the bed or Rhonegletscher in summer 2017 and 2018. Together with the spectral width of these overtones, estimations for the resonator extension can derived
Additionally to built-in python packages 'numpy', 'scipy', and 'obspy' are used.
Read in miniSEED data from crack wave events¶
The recorded pressure is saved in miniSEED format in the directory '2017events/data'. We read the raw data in and convert the digital counts to the corresponding pressure in Pa.
#read in the miniseed data for each time
import obspy
import numpy as np
p_air = 1e5 * np.exp(-2300/8435) #air pressure altitude correction at 2300 m
counts_to_bar = 1.62e5 #cnt/bar
st = obspy.read('2017events/data/crackwave*.mseed')
for trace in st:
print(trace)
trace.data = ((trace.data / counts_to_bar) - 9.39 ) * 100000 - p_air #convert in Pa and subtract air pressure
Spectral analysis of single crack waves¶
For each crack wave event we apply a similar spectral analysis. Thus we write for each step a simple function and apply it to each crack wave event that is saved in the variable st. The first function simply plots the waveform of a crack wave event.
#plot the waveform data
import matplotlib.pyplot as plt
def plot_trace(trace):
#bandpass filter the pressure data
plot_trace = trace.copy()
plot_trace.detrend('linear')
plot_trace.taper(0.01)
plot_trace.filter('lowpass', freq=20, corners=2, zerophase=False)
data = plot_trace.data
times = plot_trace.times()
#plot the waveform of a crack wave event
fig_wave, ax = plt.subplots()
plt.plot(times, data-data[0], linewidth=1)
plt.xlabel('Time (s)', fontsize=12)
plt.ylabel(r'$\Delta$ Pressure (Pa)', fontsize=12)
plt.title('Pressure oscillation during crack wave')
ax.tick_params(axis='both', labelsize=12)
fig_wave.tight_layout()
return(fig_wave)
In order to derive estimations for the resonator volume in which the crack waves travel, we have to carefully analyze the frequency content of each crack wave.
#get spectrum of crack wave event
import scipy.signal
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
def spectro(trace):
fs = trace.stats.sampling_rate
data = trace.data
#get the power spectral density
Pxx, freqs = mlab.psd(data, Fs=fs,
NFFT=min(len(data),int(fs*10)),
detrend='linear')
#get the amplitude spectral density
Pxx = np.sqrt(Pxx)
#supersample the data
freqs_super = np.linspace(0, fs/2., int(fs/2.*10))
Pxx_super = scipy.signal.resample(Pxx, int(fs/2.*10))
fig_spec, ax = plt.subplots()
plt.plot(freqs_super,Pxx_super)
plt.xscale('log')
plt.ylim([0,1.1*max(Pxx_super)])
plt.xlabel('Frequency (Hz)', fontsize=12)
plt.ylabel('Amplitude Spectral Density (Pa s)', fontsize=12)
plt.title('Spectrum of crack wave')
ax.tick_params(axis='both', labelsize=12)
fig_spec.tight_layout()
return(Pxx_super, freqs_super, fig_spec)
The crucial information for the crack wave analysis is in the spectral peaks of the Fourier spectrum. We run a peak detection function and only keep the 20 strongest peaks in the frequency spectrum.
#find peaks in spectrogram frequencies between 0 and 20 Hz
import scipy.signal
def find_peaks(Pxx, freqs):
sampl_rate = len(freqs)/max(freqs)
#only look at frequencies below 20 Hz
max_freq = 20.
max_idx = np.max(np.where(freqs <= max_freq))
freqs_small = freqs[:max_idx]
Pxx_small = Pxx[:max_idx]
#peak detection
peak_idx, _ = scipy.signal.find_peaks(Pxx_small,
height=None, threshold=None,
distance=int(0.8*sampl_rate), prominence=None,
width=0.1*sampl_rate, wlen=None, rel_height=0.5)
#only pick the 20 dominant peaks
peaks_value = Pxx_small[peak_idx]
sort_idx = np.argsort(peaks_value)
sort_idx = np.flip(sort_idx, 0)
sorted_idx = peak_idx[sort_idx]
strong_peak_idx = sorted_idx[:20]
#sort again by frequency
strong_peak_idx = np.sort(strong_peak_idx)
return(strong_peak_idx, Pxx_small, freqs_small)
In order to determine the peak position and width in the Frequency domain most reliable, we fit a gaussian to each of the maximum 20 dominant peaks.
import scipy.optimize
#define Gauss function for fitting
def gauss(x, a, x0, sigma):
return a*np.exp(-(x-x0)**2/(2*sigma**2))
#fit Gauss function to major peaks
def fit_gauss(Pxx, freqs, peak_idx):
gauss_par = np.zeros([len(peak_idx),3])
for i, idx in enumerate(peak_idx):
popt = []
p0 = [Pxx[idx], freqs[idx], 0.1] #first guess
bounds=([0.9*Pxx[idx], freqs[idx]-0.1, 0.001],[1.1*Pxx[idx], freqs[idx]+0.1, 1])
popt, _ = scipy.optimize.curve_fit(gauss, freqs[max(0,idx-5):idx+5],
Pxx[max(0,idx-5):idx+5], p0=p0, bounds=bounds)
if abs(popt[2] - 1.) > 1e-5:#append parameters only where the fit worked
gauss_par[i,:] = popt
gauss_par = gauss_par[~np.all(gauss_par == 0, axis=1)]
return(gauss_par)#amplitude, peak position, width
With a simple plot we can check the reliability of the gauss fits.
#plot spectrum of crack wave with the fitted gauss functions
import matplotlib.pyplot as plt
def plot_gauss(Pxx, freqs, peaks_final):
fig_gauss, ax = plt.subplots()
plt.plot(freqs, Pxx, linewidth=2, color='black')
for fit in peaks_final:
plt.plot(freqs, gauss(freqs,*fit))
plt.xlim([0,20])
plt.ylim([0,1.1*max(Pxx)])
plt.xlabel('Frequency (Hz)', fontsize=12)
plt.ylabel('Amplitude Spectral Density (Pa s)', fontsize=12)
plt.title('Spectrum with fitted gauss functions')
ax.tick_params(axis='both', labelsize=12)
fig_gauss.tight_layout()
return(fig_gauss)
With the code below, we run the above functions for each crack wave event. We save the plots of the pressure waveform, spectra and gauss fits together with the csv files including the fitted peak parameters separately for each event. As an example we show the plots of the first crack wave event below.
import os
from obspy.core import UTCDateTime
for i, trace in enumerate(st):
time_str = UTCDateTime.strftime(trace.stats.starttime, format='%Y-%m-%d-%H-%M-%S')
#output of the spectral analysis will be saved here:
savepath = '2017events/analysis/spectra/'+time_str+'/'
if not os.path.exists(savepath):
os.makedirs(savepath)
#plot waveform
fig_wave = plot_trace(trace)
fig_wave.savefig(savepath+time_str+' wave.pdf')
#calculate and plot amplitude spectral density (ASD)
Pxx, freqs, fig_spec = spectro(trace)
fig_spec.savefig(savepath+time_str+' spectro.pdf')
#find peaks in the ASD
strong_peak_idx, Pxx_small, freqs_small = find_peaks(Pxx, freqs)
#fit a gaussian function to the peaks and check the fit
gauss_par = fit_gauss(Pxx_small, freqs_small, strong_peak_idx)
#peaks_final = check_fit(gauss_par)
#plot the fits
fig_gauss = plot_gauss(Pxx_small, freqs_small, gauss_par)#peaks_final)
fig_gauss.savefig(savepath+time_str+' gauss.pdf')
#reorder and save frequency peak parameters to csv file
peaks_final_reorder = gauss_par[:,[1,2,0]] #Frequency, Width, Amplitude
peaks_final_reorder[:,1] *= 1.67 #convert sigma in width at 1/sqrt(2)
np.savetxt(savepath+time_str+'_peaks.csv',
peaks_final_reorder, delimiter=", ", fmt=['%.2f','%.2f','%.2f'] ,
header='Frequency (Hz), Width (Hz), Amplitude (Pa s)')
if i == 0:
plt.show()
plt.close('all')
Kernel Density Estimation of Peaks¶
After we determined the spectral peak positions and widths for each crack wave, we use this complete ensemble to calculate most probable peak positions and widths. To do so, we first read in all the peak positions, widths and amplitudes from each crack wave event.
#output from the kernel density estimation will be saved here:
import os
savepath = '2017events/analysis/kde/'
if not os.path.exists(savepath):
os.makedirs(savepath)
#read in the peak csv files and create one list with all frequencies, widths and amplitudes
import glob
import csv
import numpy as np
files = glob.glob('2017events/analysis/spectra/*/*_peaks.csv')
freq_list = np.array([])
width_list = np.array([])
amp_list = np.array([])
for file in files:
with open(file, newline='') as csvfile:
readCSV = csv.reader(csvfile, delimiter=',')
next(readCSV, None)
for row in readCSV:
freq = float(row[0])
width = float(row[1])
amp = float(row[2])
freq_list = np.append(freq_list, freq)
width_list = np.append(width_list, width)
amp_list = np.append(amp_list, amp)
sort_idx = np.argsort(freq_list)
freq_list = freq_list[sort_idx]
width_list = width_list[sort_idx]
amp_list = amp_list[sort_idx]
Now we apply a kernel density estimation (KDE) with a bandwidth of 0.3 Hz on the the entirety of all spectral peak positions.
#get the kernel density estimation (KDE) from the peak distribution
from sklearn.neighbors import KernelDensity
import numpy as np
freqs = np.linspace(0, 20, 200)
#fit kernel density model on the data
kernel = 'epanechnikov'
kde = KernelDensity(kernel=kernel, bandwidth=0.3).fit(freq_list[:, np.newaxis])
#evaluate density model on the data, outputs array of log(density)
log_dens = kde.score_samples(freqs[:, np.newaxis])
kde_appl = np.exp(log_dens)
We plot the calculated KDE. It essentially shows the density of the spectral peak distribution. The small black tickes along y=0 are the individual peak positions.
#plot the KDE
import matplotlib.pyplot as plt
fig_kde, ax = plt.subplots()
ax.plot(freqs, kde_appl)
ax.plot(freq_list, np.zeros((len(freq_list))), '+k')
ax.set_ylim(-0.01, 1.1*max(kde_appl))
plt.xlabel('Frequency (Hz)', fontsize=12)
plt.ylabel('Probability density', fontsize=12)
plt.title('KDE of peak positions')
ax.tick_params(axis='both', labelsize=12)
fig_kde.tight_layout()
fig_kde.savefig(savepath+'kde.pdf')
plt.show()
Similar to the above spectral analysis we derive the peak positions and corresponding widths of the KDE by fitting gaussians. We save the output into a csv file and plot the fits into KDE.
#get peaks of the KDE
kde_peak_idx, kde_appl, freqs = find_peaks(kde_appl, freqs)
#fit gauss to peaks
kde_gauss_par = fit_gauss(kde_appl, freqs, kde_peak_idx)
#reorder the peaks and save it to a csv file
kde_gauss_par_reorder = kde_gauss_par[:,[1,2,0]] #reorder as Frequency, Width, Amplitude
np.savetxt(savepath+'kde_peaks.csv',
kde_gauss_par_reorder, delimiter=", ", fmt=['%.2f','%.2f','%.2f'],
header='Frequency (Hz), Sigma (Hz), Probability Density')
#plot the gauss fit
fig_kde_gauss, ax = plt.subplots()
plt.plot(freqs, kde_appl, color='black')
for fit in kde_gauss_par:
plt.plot(freqs, gauss(freqs,*fit))
plt.xlim([0,20])
plt.ylim([0.001,1.5*max(kde_appl)])
plt.xlabel('Frequency (Hz)', fontsize=12)
plt.ylabel('Probability density', fontsize=12)
plt.title('KDE of peak positions with fits of gaussians')
ax.tick_params(axis='both', labelsize=12)
fig_kde_gauss.tight_layout()
fig_kde_gauss.savefig(savepath+'kde_gauss.pdf')
plt.show()
From the peak position and peak widths of the KDE of all spectral peak positions, we select those spectral peaks that lie within one sigma of the peak positions determined in the KDE. From this subset, we calculate the mean of their spectral peak positions and their spectral widths as well as the standard deviation.
#extract within one sigma all spectral peak frequencies and widths
#and calculate their mean and standard deviation
arr = np.zeros([len(kde_gauss_par_reorder),5])
for i, (frequency, sigma, kde) in enumerate(kde_gauss_par_reorder):
fn_idx = [(freq_list > (frequency-sigma)) & (freq_list < (frequency+sigma))]
fn_freq = freq_list[fn_idx]
fn_width = width_list[fn_idx]
fn_freq_mean = np.mean(fn_freq)
fn_freq_std = np.std(fn_freq)
fn_width_mean = np.mean(fn_width)
fn_width_std = np.std(fn_width)
n_peaks = len(fn_freq)
arr[i] = [n_peaks, fn_freq_mean, fn_freq_std, fn_width_mean, fn_width_std]
#save to csv file
savepath = '2017events/analysis/final/'
if not os.path.exists(savepath):
os.makedirs(savepath)
np.savetxt(savepath+'final_peaks.csv',
arr, delimiter=", ", fmt=['%.0f','%.2f','%.2f','%.2f','%.2f'],
header='Number of Peaks, Frequency Mean, Frequency Std, Width Mean, Width Std')
Calculation of physical crack parameters¶
With the derived spectral peak positions and their widhts in the frequency domain, the dimensions of the basal water sheet patch where we recorded crack waves, can be estimated.
#calculate the crack length and aperture from the frequency peaks and widths
#only take the first 5 frequency peaks into account
fs = arr[:5,1] #frequencies
fs_err = arr[:5,2] #frequency errors
dfs = arr[:5,3] #widths
df_err = arr[:5,4] #width errors
nu = 1.7e-6 #m^2/s kinematic viscosity of water
rho = 1000 #kg/m^3 density of water
pi = 3.14
#do the calculation for each of the frequency peaks
all_vals = np.zeros([len(fs),7])
for i, (f,ferr,df,dferr) in enumerate(zip(fs, fs_err, dfs, df_err)):
Ls = ([])
ws = ([])
Lserr = ([])
wserr = ([])
Q = f / df #Quality factor describing the attenuation of the signal
Q_err = Q * np.sqrt((ferr/f)**2 + (dferr/df)**2)
for icerock in ['ice', 'rock']: #['ice'], ['rock'] for only ice and only rock
if icerock == 'ice':
nu_s = 0.35 #poisson's ratio for ice
G = 3.8e9 #Pa shear modulus for ice
elif icerock == 'rock':
nu_s = 0.25 #poisson's ratio for rock
G = 27e9 #Pa shear modulus for rock
G_star = G/(1-nu_s) #elastic plain strain modulus
L = 0.5 * (pi * nu * (G_star/rho)**2 * Q**2/f**5)**(1/6) #m crack length
w = Q * np.sqrt(nu/(pi * f)) * 1000 #mm full crack aperture
Ls.append(L)
ws.append(w)
L_err = L * np.sqrt(0.5**2 * (ferr/f)**2 + (1/3)**2 * (dferr/df)**2) #m crack length error
w_err = w * np.sqrt(0.5**2 * (ferr/f)**2 + (dferr/df)**2) #mm crack aperture error
Lserr.append(L_err)
wserr.append(w_err)
L_mean = np.mean(Ls)
L_errmean = np.mean(Lserr)
w_mean = np.mean(ws)
w_errmean = np.mean(wserr)
all_vals[i] = [f, Q, Q_err, L_mean, L_errmean, w_mean, w_errmean]
np.savetxt(savepath+'final_values.csv',
all_vals, delimiter=", ", fmt=['%.1f','%.1f','%.1f','%.1f','%.1f','%.1f','%.1f'],
header='Peak Frequency, Q, Q Std, Length (m), Length Std (m), Aperture (mm), Aperture Std (mm)')
print('Extension of the basal water sheet patch:\n'
'Length: {0:.1f} m'.format(all_vals[0,3]),'\n'
'Length Std: {0:.1f} m'.format(all_vals[0,4]),'\n'
'Aperture: {0:.1f} m'.format(all_vals[0,5]),'\n'
'Aperture Std: {0:.1f} m'.format(all_vals[0,6]))