Repository associated with the manuscript "On the quieting of glaciohydraulic tremor" by Chmiel et al.

Data products and code for: "On the quieting of glaciohydraulic tremor" Authors: M. Chmiel, N. Caldera, F. Gimbert, G. Olivier, A. Guadagnini, D. Gräff, M. Köpfli, F. Walter Submitted to: Earth and Planetary Science Letters Date: June 14, 2025 Overview This repository contains the data products, MATLAB functions, and Jupyter Notebooks required to reproduce the results and figures from the manuscript "On the quieting of glaciohydraulic tremor" by Chmiel et al. Repository Contents Code Bartlett_xyz_grid.m This function implements a Bartlett beamformer to map potential source locations across a predefined 3D grid. It processes data from a sensor array at various frequencies to create a complete spatial power map. Bartlett_xyz_min.m This function is designed for use with an optimization algorithm (like MATLAB's fminsearchbnd) to efficiently find the location of the single dominant source. It calculates the Bartlett beamformer power for a single 3D point (xyz) and returns its negative value. An optimization routine minimizes this function's output, which in turn maximizes the beamformer power. Compute_hydraulic_potential.m This script loads Digital Elevation Model (DEM) data for the Rhone Glacier, calculates the hydraulic potential, and visualizes the results. MFD_UpstreamArea.m This script computes the upstream area as the sum of all hydraulically connected upstream grid cell areas using the multiple-direction flow accumulation method. Model_channel_discharge_as_a_function_of_filling_depth.ipynb This Jupyter Notebook calculates changes in conduit discharge (Q) and hydraulic radius (R) as a function of conduit filling for different subglacial conduit geometries. Data products MFP_results /01_gridsearch/ This folder contains the output from our unconstrained 3D grid searches (example Figure 3B in the manuscript). These results were used to confirm that tremor sources originate near the glacier bed, justifying the use of a depth-constrained search. Each directory, 10-20 and 6-8 (for the corresponding frequency bands in Hz), contains daily folders with Matched Field Processing (MFP) results calculated over 5-minute intervals. Each *.mat file contains a structure named “Localisations” with the following information: vfinal: 1800 [m/s] final optimized phase wave velocity xx: Easting coordinates for the grid search yy: Northing coordinates for the grid search zz: Altitude coordinates for the grid search vect_data_adaptive_per: [8×20×30 double] - complex vector data that contains phase information (dimension: N_stations x N_frequencies x N_segments) freq_in: A vector containing the frequency vector for which MFP is calculated BA_value: Maximum Bartlett output value obtained through optimization xyz_opt: XYZ coordinates of the optimal source location through an optimization approach (see also /02_depth_constrained_search/). /02_depth_constrained_search/ This folder contains the final, quality-controlled results used for interpretation and to generate Figures 1B-C and 5A in the manuscript. 10_20_MFP_locations.csv: The complete catalog of source locations with MFP from the depth-constrained analysis in the 10-20 Hz band. This file includes quality flags. XYZ-coordinates of the starting point for minimization = [2672320, 1161360, 2314] 6_8_MFP_locations_four_starting_points.csv: The complete catalog of source locations with MFP from the depth-constrained analysis in the 6-8 Hz band. This file includes quality flags. XYZ-coordinates of the four starting points for minimization: [2672250, 1161260, 2314; 2672250, 1161380, 2314; 2672370, 1161260, 2314; 2672370, 1161380, 2314] Seismic stations used: 3 surface stations (RA81, RA82, and RA87) and 5 borehole stations (RA91A, RA91B, RA92A, RA93A, and RA93B). /03_verification_full_array_search/ This folder contains results from the unconstrained verification run using the full seismic array with multiple starting points for the 6-8 Hz frequency band, as shown in Figure E.6. As described in the manuscript, this run confirms a near-bed source origin, though with lower horizontal precision than our primary results. 6_8_MFP_locations_16_starting_points_whole_array.csv: 3D locations from the verification runs using the entire array: surface stations RA81-RA88 and deep-borehole stations RA91-RA93. Coordinates of the starting points for minimization: [2672150, 1161190, 2314; 2672150, 1161310, 2314; 2672150, 1161430, 2314; 2672150, 1161550, 2314; 2672270, 1161190, 2314; 2672270, 1161310, 2314; 2672270, 1161430, 2314; 2672270, 1161550, 2314; 2672390, 1161190, 2314; 2672390, 1161310, 2314; 2672390, 1161430, 2314; 2672390, 1161550, 2314; 2672510, 1161190, 2314; 2672510, 1161310, 2314; 2672510, 1161430, 2314; 2672510, 1161550, 2314] Xcorr_selection This section contains examples of cross-correlation functions between borehole-borehole stations (station pair 91A-93A, component ZZ shown in Figures 3A and E.5) and borehole-surface stations (station pair RA82-RA91, component ZZ). The MATLAB files contain monthly cross-correlation functions. Each file includes the following data: stack_wav: Monthly cross-correlation functions sampled at 500 Hz, with dimensions 8640x8001 (N_correlations x Time samples in the time lag vector). The correlations start on July 22, 2020, and each correlation function corresponds to a 5-minute long time window. time_lag: A time lag vector between -8s and 8s, with a time step of 0.002 s. Dispersion_curves_slowness_2modes.txt: This file contains dispersion curves of the fundamental mode (Mode 0) and the second higher mode (Mode 2) in the format: Frequency (Hz), Slowness (s/m). The corresponding dispersion measurements are plotted in Figure 5A.