Diurnal expansion and contraction of englacial fracture networks revealed by seismic shear wave splitting

Fractures contribute to bulk elastic anisotropy of many materials in the Earth. This includes glaciers and ice sheets, whose fracture state controls the routing of water to the base and thus large-scale ice flow. Here we use anisotropy-induced shear wave splitting to characterize ice structure and probe subsurface water drainage beneath a seismometer network on an Alpine glacier. Shear wave splitting observations reveal diurnal variations in S-wave anisotropy up to 3%. Our modelling shows that when elevated by surface melt, subglacial water pressures induce englacial hydrofractures whose volume amounts to 1-2 percent of the probed ice mass. While subglacial water pressures decrease, these fractures close and no fracture-induced anisotropy variations are observed in the absence of meltwater. Consequently, fracture networks, which are known to dominate englacial water drainage, are highly dynamic and change their volumes by 90-180 % over subdaily time scales. Daily variations in seismic S-wave anisotropy reveal the dynamic response of englacial fractures to changing meltwater supply, according to shear wave splitting measurements from a seismometer network on the Rhonegletscher, Switzerland.

L iquid water controls the fast flow of glaciers and ice sheets. Meltwater feeds the subglacial hydraulic system, where it modifies pressures to reduce or enhance ice sliding 1 . The englacial drainage system plays an important role in glacial hydraulics and consequently water-controlled basal sliding. It provides water storage space, hydraulically links the ice surface to the bed and contains different forms of water flow pathways 2 . Of these, networks of linked fractures are believed to dominate over efficient ice-walled conduits 3 . Englacial fracture networks also facilitate glacier surging 4 and allow for catastrophic drainages of supraglacial 5 and ice-marginal 6 lakes. These drainage events increase ice sheet flow on day-long 5 and multiannual 7 time scales.
Our knowledge of englacial drainage relies primarily on pointlike borehole observations 3,8 . Radar echo sounding and active seismic surveys also offer images of englacial and subglacial water inclusions 9,10 . However, required labour limits repeated surveys to few time snapshots 10 or confines them to isolated locations on the ice 11,12 . This complicates monitoring the englacial drainage system at seasonal and daily time scales over which hydraulic states and water pressures within glacier ice can change significantly 13 . Accordingly, there exists little information about the temporal variability of englacial fracture networks beyond indirect inferences from lake drainage events 14 .
Here we analyse shear wave splitting (SWS) in seismic signals from icequakes to monitor changes in the englacial drainage system of Rhonegletscher, an Alpine glacier in Switzerland (Fig. 1). We model the effects of preferentially aligned ice crystals 15,16 and preferentially aligned englacial fractures 17 to identify changes in water-filled fracture volume. Applied to naturally occurring seismicity at the glacier base, this approach describes the englacial fracture state at an unrivalled temporal resolution.
Results and discussion SWS in temperate glacier ice. We rely on seismic SWS 18 as a direct indicator of elastic anisotropy. SWS occurs when a shear wave enters an anisotropic medium splitting the shear wave into two orthogonally polarized phases propagating at different velocities. The underlying anisotropy may result from preferential alignment of crystals 15 ("crystal orientation fabric", COF) or fractures 18,19 within the propagation medium. For COF, the elastic properties of crystals are directionally dependent. For fracture-induced anisotropy, shear wave components whose polarization axis does not cross fractures propagate faster than shear wave components with fracture-crossing polarization 20 (Fig. 2a).
SWS has been exploited in numerous settings to measure the anisotropy of the Earth's subsurface 18 ranging from laboratory 21 to crustal and mantle scales 19 . Highly sensitive to fracture detection, SWS is an effective tool in exploration contexts providing reservoir characterization in geothermal 20 as well as oil and gas 22 applications. SWS constitutes a means for monitoring the density, geometry, and connectivity of fractures in hydraulically stimulated shale gas reservoirs 23,24 .
We investigate temporal changes in anisotropy of glacier ice revealed by SWS. Our 2018 study site is located in a high-melt zone of Switzerland's Rhonegletscher (Fig. 1a), one of the world's Locations of basal stick-slip icequake clusters (red stars) and seismometers (blue triangles) active/operational during summer used in this study. Grey coloured items indicate the positions of remaining summer clusters/ seismometers that were discarded due to low-quality data or poor temporal coverage. Ruby pins indicate differential GPS stations. Circled blue triangles indicate seismometers that were operational also during winter. Circled red star denotes locations of three nearly co-located through different stick-slip clusters: one active only during winter and two active during summer. The dashed black line approximates the central flow line. c Thin section of the ice core (images obtained from cross-polarized light microscopy, the location shown in (a)) with different colours depicting crystals with different c-axis orientations. The two examples at 45 m depth show fissures with small crystals (shown in dark blue) grown during refreezing upon core extraction. Red dotted lines indicate the mutual orientation of vertical (rectangular) and horizontal (circular) cuts. The ice flow direction is shown with black arrows. In both examples, the orientation of the fissures is 45°with respect to glacier flow and thus aligned with the nearby crevasses (marked with black dashed lines). Source of a, b base maps: Swiss Federal Office of Topography. most extensively studied temperate glaciers 10,25 . At the study site, the ice thickness reaches 200 m as determined from interpolation of radar transects 26 and hot-water drilling 27 and the glacier exhibits microseismic stick-slip sliding over its bed 28 . We analyse stick-slip icequakes recorded with arrays of shallow (1-3 m in-ice depth) borehole seismometers (LE-3D/BHs by Lennartz, flat frequency response between 1 and 100 Hz) installed above stickslip asperities (Fig. 1b). These icequakes tend to cluster at asperities where frictional properties of the ice-bed interface lead to elastic strain accumulation, which is subsequently released via seismogenic slip 28,29 .
Stick-slip icequakes were recorded within one winter month (13.02-20.03.2018) and one summer month (21.07-22.08.2018). Summer data were recorded using nine seismometers, of which three were operational also during winter supplemented by four additional geophones. Signal detection was based on typical waveform characteristics and cross-correlation searches and stick-slip sources were located with a probabilistic travel time inversion (Methods and Gräff and Walter 29 ). We use winter data as a reference and focus on analysing summer data among which we identify 1158 basal icequakes, whose signals have nearly identical waveforms 28 , distributed among seven clusters. A typical stick-slip seismogram lasts around 0.1 s and is dominated by Pand S-waves with frequencies up to a few hundred Hz (Fig. 2b).
We focus on the shear waves of our stick-slip icequakes, which englacial anisotropy is expected to split into two orthogonally polarized phases travelling with different velocities 30 (Fig. 2a). From measurements of the polarization direction of the fast S-wave phase and the delay time between the fast and slow phases (Fig. 2c, d), both the fast axis orientation of the medium and the magnitude of anisotropy along a given ray path can be reconstructed 19 .
SWS measurements were taken for each event at each basal stick-slip cluster-seismometer pair using the automated technique by Wuestefeld 31 . This processing resulted in a catalogue of delay times and fast shear wave polarization angles together with associated uncertainties based on the maximum differences of each parameter values among different analysis time windows.
Only high-quality SWS measurements were accepted, which limited the analysed stick-slip summer data to three clusters. The findings presented in this study are contingent upon this quality control ( Supplementary Fig. S1) whose criteria include: high signal-to-noise ratio, effective ray-frame rotation, a single minimum of error surface, consistent result regardless of the analysis window length, and low measurement uncertainties (Methods). We analyse the 11-day long time period covered quasi-continuously by 472 measurements. Delay times were normalized by source-receiver distance and mean S-wave velocity to obtain the magnitude of S-wave anisotropy in percent (Methods).
Time series of S-wave anisotropy magnitude obtained for winter and summer datasets differ substantially. During winter, individual icequake cluster-seismometer pairs show constant magnitudes of S-wave anisotropy over time (Fig. 3a) and these magnitudes depend on the orientation of the ray path in the anisotropic medium. In contrast, during summer, magnitudes of S-wave anisotropy vary diurnally and correlate with meltwater discharge (Fig. 3b) measured 3 km downstream from the terminus of Rhonegletscher. Measurement details for individual icequake cluster-seismometer pairs active in both winter and summer seasons (a subset of data in Fig. 3) are given by Supplementary Fig. S2. The fast S-wave polarization angles also exhibit daily variations during the summer period, only ( Supplementary Fig. S3).
The spatial clustering and nearly identical waveforms of our basal stick-slip icequakes are ideal for monitoring temporal changes in the medium traversed by the seismic rays. SWS depends solely on the anisotropic features of the medium aggregated over the whole ray path, is independent of seismic source mechanisms or medium changes outside the ray paths 18 , and is not subject to free-surface effect ( Supplementary Fig. S4). Therefore, the presented SWS reflects effective anisotropy of glacial ice from the bed to the surface. For our stick-slip icequakes, reflected and refracted signals may traverse a thin basal till layer 32 and possibly contribute to observed SWS due to variations in elastic till properties resulting from pore pressure fluctuations. However, we exclude possible till contribution to SWS because low till seismic velocities and seismic network geometry imply near-vertical rays limiting horizontal propagation of shear waves through the till lenses. Thus, we attribute the observed diurnal variations in S-wave anisotropy to periodic and therefore reversible changes within the ice.
Source of anisotropy variation. Crystal orientation fabric (COF) is known to cause elastic anisotropy in glacial ice 15,16 but does not explain observed diurnal variations in S-wave anisotropy. Change in COF results from viscoelastic deformation and dynamic recrystallization of the ice governed by the local stress regime 33 . However, dynamic recrystallization needed to substantially change fabric occurs after some 10% of strain 34 which at stimulated strain rates in laboratory simulations involves at least days 35 . Although COF may contribute to seismic anisotropy, it does not explain the observed diurnal changes.
Our observed anisotropy variations are best explained in terms of changes in preferentially aligned englacial fractures such Fig. 2 Shear wave splitting mechanism with an example from icequake data. a Principle of fracture-induced shear wave splitting (SWS), adopted after Rial et al. 20 . b An example of a Rhonegletscher stick-slip icequake seismogram in ray-frame coordinates, the highlighted area indicates the S-phase presented in (c, d). c Time delay between two S-wave phases before and after SWS correction. d Relative amplitude of pre and post-correction particle motion of signals shown in (c).
as small millimetre to centimetre-scale cracks or larger, meterscale crevasses. Pressurized water within the glacier can induce fractures ("hydrofractures") 36 and pressure variations enhance opening 37 and closing 38 of such fractures. Englacial fractures tend to align in preferred orientations relative to principal strain axes. Such fracture patterns slow down S-waves, whose particle motion crosses the fractures 17 . Within the efficient subglacial drainage system, water pressure variations beneath high-melt zones reach hundreds of metres of water column and during peak pressure the overlying ice approaches flotation 13 . Melt-induced pressure variations thus explain our observed S-wave anisotropy in terms of preferentially oriented hydrofractures, which is analogous to seismic anisotropy changes in geothermal reservoirs during hydraulic stimulation 20 .
Englacial fractures connected to the efficient subglacial drainage system are filled and/or pressurized with water during warm daytimes with increased surface runoff. Consequently, in dry or semi-filled spaces the air is replaced with water, which may affect the anisotropy, both magnitude and polarization 23,39 . However, our scenario requires fractures to be efficiently linked with water pathways allowing the fluid pressure to equilibrate while the seismic wave is passing (low-frequency condition). That makes the fluid substitution effect less influential than in the case of isolated fractures or highly viscous fluids 40 . Our numerical modelling (Methods) shows that the expected change in the S-wave anisotropy magnitude due to air-to-water fluid substitution for ray paths available in field data span from −0.15% to 0.1% ( Supplementary Fig. S5), and hence, is an order of magnitude weaker than the influence of cracks. This change may be further reduced since our study site is located in a highmelt zone meaning that cracks will unlikely drain completely. Consequently, mere changes of water fill state, pressure or rearrangement within englacial fractures do not explain our observed anisotropy variations.
Diurnal temperature changes at the surface can penetrate the ice via heat conduction but attenuate to 5% at a depth of only 50 cm 41 . In contrast, diurnal changes in englacial water content driven by the day/night surface melt cycle can affect large parts of the ice column 12 . The maximum anisotropy during high-melt daytimes and the nighttime minimum (Fig. 3b) are in agreement with hydrofracture formation during the day and fracture network closure during the night. Compared to hydraulically opened fractures, closed fractures transfer stresses across the touching fracture surfaces, which reduces displacement discontinuities and thus the fracture-related component of the elastic compliance tensor 42  Modelling fracture-induced anisotropy. We calculate density and strike (azimuth of a horizontal axis within the fracture plane) of fracture networks with preferred orientation during low-and high-anisotropy time periods (Fig. 3b). The low magnitude of anisotropy during nighttimes characterizes a "background anisotropy" comprised of some combination of ice COF and fractures remaining open even during low water pressures. We select 8 h-long time windows that cover low anisotropy measurements (grey boxes in Fig. 3b).
For fractures much smaller than S-wavelengths (19 m at a frequency of 100 Hz and a velocity of 1900 m s −1 ), anisotropy does not depend on fracture size. This assumption is met for millimetre to centimetre scale cracks but violated for crevasses with meter-wide opening and lateral extents of tens of metres. To model seismic anisotropy induced by millimetre to centimetre scale cracks, we invoke elastic effective medium theory 43 . We extend the approaches of Harland et al. 44 and Smith et al. 15 to invert SWS data (delay times and fast polarization angles) for ice fabric (combinations of COF) and set of aligned vertical pennyshaped cracks 45 saturated with water 46 . To this end, we use the additional compliance approach of Schoenberg and Sayers 47 as was done by Verdon et al. 48 in the SWS context. We consider different combinations and proportions of any two ice fabrics including: cluster, thick girdle, partial girdle, and horizontal partial girdle with different opening angles and orientations (defined as in Maurel et al. 49 ) superimposed with cracks parametrised by density and strike. In agreement with borehole observations of predominant steeply dipping fractures 3,8 including nearby ice core samples (Fig. 1c), the superimposed cracks have a vertical orientation indicating nearly horizontal principal stresses.
The best-fitting low-pressure background model is composed of a mixture of 80% of thick girdle fabric with an opening angle of θ = 60°oriented at azimuth 130°CWN (clockwise from North) and 20% of the partial girdle with a narrow opening angle of θ = 10°oriented at 180°CWN with 0.9% of vertical cracks (% of total ice volume) striking at 180°CWN (Fig. 4a). A second model run determines the best fit between peak pressure anisotropy data (yellow boxes in Fig. 3b) and previously inverted background ice fabric supplemented with another set of vertical cracks. In addition to the low-pressure background model, the best fitting model contains another 0.8% of cracks striking at 160°CWN (Fig. 4b), which constitutes a 90% rise in crack density. We repeat the peak pressure anisotropy inversion, but this time excluding three days when meltwater discharge was relatively low and did not exceed 10 m 3 s −1 (hatched yellow boxes in Fig. 3b). During these days, hydrofracturing is expected to be reduced. With respect to the low-pressure background model, this results in 1.6% of additional cracks striking at 180°CWN (Fig. 4c), which constitutes a 180% rise in crack density in respect to the low anisotropy scenario. Newly developed fractures are aligned with the fast direction of the background model (Fig. 4). This is expected since daily variations in the magnitude of S-wave anisotropy (Fig. 3) are larger than variations in fast polarization angles ( Supplementary Fig. S3).
The fracture volume changes occur over diurnal time scales much shorter than the 11 days over which previous borehole observations documented fracture opening 8 . Fracture networks thus have to be well connected to the efficient subglacial drainage system as fracture opening and closure occurs during meltinduced high and low-pressure periods, respectively.
The presented low-pressure anisotropy model has 9 parameters and is thus less constrained and more prone to parameter tradeoff than the 2-parameter (only density and strike of extra fractures) peak-pressure models. Also, the COF of a temperate alpine glacier may be more complex and spatially variable. Hence, our model may not capture the complete ice fabric description. Nevertheless, it serves its purpose as a reference to track fracture intensification during high anisotropy periods since we associate all the increase in anisotropy with fracture growth. The modelled storage volumes still are in agreement with previous measurements of water-filled englacial void spaces with volumes of 0.5-3% of the glacier region 8 . This supports our finding that the additional 0.8-1.6% of cracks corresponds to a ca. 90-180% increase in storage volume with respect to the 0.9% crack density during low-pressure times.
During both high-and low-melt periods, the inverted crack strike is nearly parallel to the ice flow direction (Fig. 4). Lack of crevasses and decreasing distances between GPS stations during most of the summer period indicate horizontal compression at the study site ( Supplementary Fig. S6) Diurnal variations in GPS stations separation suggest that diurnal surface strain variations in both longitudinal and transverse directions are more than two orders of magnitude below the expected volume changes of 0.8-1.6% (Supplementary Fig. S7). Furthermore, the ice surface shows no major crevasses and linear features visible in airborne imagery are likely healed relict crevasses, which were opened during advection through steep terrain upglacier ( Supplementary  Fig. S8). These observations show that the dynamic fracture network does not break the ice surface and consequently, the 90-180% fracture volume variations, which represent an average over steep ray paths, are exceeded at some depths.
Water-routing fractures, which terminate a few metres below the ice surface have been found in other glacial environments 50 . In our case, the water, which occupies the additional 90-180% fracture volume during peak melt, is too large to be produced locally, given maximum daily surface ablation of 10 cm (Landmann et al. 51 ) equivalent to only~0.05% of local ice volume. Therefore, the fracture network must store additional water from surrounding regions, potentially in relation to local bed overdeepening 28 .
We expect millimetre to centimetre-scale cracks to be responsible for the changing S-wave anisotropy. Meter-scale scale crevasses are also influenced by hydrofracturing at their tips 36 , but they remain open even when convected into compressional flow regimes, which testifies to longer time scales required for crevasse closure 52 . On different glaciers, borehole studies found typical englacial fracture opening of a few centimetres 3,8 , and smaller cracks likely evade borehole observations. In ice cores on Rhonegletscher's terminus, vertical cracks of a few centimetre length and about half a centimetre width were found at various depths ( Fig. 1a and c).
Frequently occurring repetitive natural seismic sources were leveraged to resolve the temporal fracture-driven variations in englacial anisotropy. Despite relatively inefficient drainage, Fig. 4 Shear-wave splitting signature resulting from different ice-fabric models. a-c Upper hemisphere projections of measured SWS data (heavy coloured ticks with white outlines) and the best fitting ice-fabric models (black ticks and background colour). The tick position indicates the azimuth and inclination (radial coordinate) of the ray path, the tick orientation indicates the fast S-wave polarization, colours and lengths of tick marks indicate the magnitude of S-wave anisotropy, and a thick black arrow indicates strike of vertical fracture set. a Inversion result of the low anisotropy scenario (SWS data within grey boxes in Fig. 3). b Inversion results for high anisotropy scenarios with all available SWS data (within all yellow boxes in Fig. 3) resulting in 0.8% of extra fractures striking at 160°. c Same as in (b) but using only days with elevated meltwater discharge (more than 10 m 3 s −1 ; data within not hatched yellow boxes in Fig. 3)  fracture networks transport more water than other components of the englacial drainage system 3 . Our study shows that the volume of this dominant water pathway is highly dynamic and reacts to changes in the surface melt, which modify subglacial water pressures.
The identified changing fracture volumes indicate that hydrological models must not consider englacial storage as a constant unit. Storage volume instead depends on subglacial water pressures and thus on seasonal variations in the configuration of the subglacial drainage system. Similarly, one-time anisotropy measurements cannot be assumed to capture fabric needed to describe ice rheology in numerical ice sheet models but include an unknown fracture component. Hydrofracturing thus changes the bulk englacial structure on short time scales that have yet to be considered in glacier hydraulic and dynamic research.

Methods
Icequakes detection and location. We use three-component Lennartz 3D/BHs seismometers drilled 1-3 m into the ice and deployed close to the central flow line of the glacier. The sensor's eigenfrequency is 1 Hz; the response is flat up to 100 Hz and the sampling frequency is 500 Hz. For the basal icequakes detection, we search for high-frequency content signals in the seismograms and check for pure P-and S-wave existence (Fig. 2b). On the found icequakes we run a hierarchical cluster analysis based on cross-correlation similarity values. Within the formed clusters, we stack the single icequakes and run a cross-correlation search on the continuous data with the stacks for each cluster. We pick the first brakes of the direct P-and S-waves with an uncertainty of one sample. To determine hypocentres, we use a probabilistic non-linear hypocentre earthquake location search approach 53 that accounts for picking and velocity model uncertainties. We use a homogenous isotropic velocity model with V P = 3750 (±5%) m s −1 and V S = 1875 (±5%) m s −1 . Details of the detection and location procedure are described in Gräff and Walter 29 .
SWS measurement methodology. SWS data processing aims at determining two parameters: time delay dt which is proportional to the magnitude of the anisotropy and fast polarization angle φ determined by the symmetry of the anisotropic medium. Recorded waveforms are rotated to the ray-frame coordinates to minimize P-wave energy on the S 1 and S 2 components (Fig. 2b). Next, the SWS correction that removes the effect of anisotropy (the splitting) is applied by rotating the two S components by the fast polarization angle φ and shifting their relative arrival time by dt. Consequently, two matching waveforms should be found (Fig. 2c) resulting in the linearization of the particle motion in a fast-and slow-axis coordinate system (Fig. 2d). The SWS correction parameters are determined by performing a grid search over all possible fast S-wave polarizations and time delays and then retrieving the minimum error solution from the resulting error surface. The best solution can be found using either the second eigenvalue of the covariance matrix in a predefined shear-wave window 54 or waveform cross-correlation 55 . As the result depends on the position and length of the analysis window, the operation is performed for various-length S-wave windows to provide a stable solution regardless of window choice 56 . Uncertainties for fast polarization angle and delay time are determined based on the maximum differences of each parameter value among different analysis time-windows. To deal with a large number of measurements, we use the automated method proposed by Wuestefeld et al 31 . that utilizes both eigenvalue and cross-correlation methods. Filtering and error analysis provide further steps to identify unique results as typically only 10-20% of the source-receiver combinations in a given dataset produce reliable SWS measurements, with noise and unfavourable ray directions usually responsible for unreliable measurements 57 . A quality assessment is performed using a machine learning algorithm trained for our dataset using a subset of measurements evaluated manually following acceptance criteria defined by Teanby et al. 56 . These criteria include: high signal-to-noise ratio, effective ray-frame rotation, a single minimum at error surface, consistent result regardless of the analysis window length, and low measurement uncertainties (visualized in Gajek et al. 24 ). Low-quality measurements are discarded.
Lastly, the magnitude of S-wave anisotropy A (%) is calculated using: where V Smean is the mean S-wave velocity (1875 m s −1 ), dt is measured delay time, and r is the distance between source and receiver.
Crack-filling fluid substitution. We modelled fluid substitution from filled with air (dry) to water-saturated conditions for ice by applying the Brown-Korringa fluid substitution formula to dry ice with preferentially oriented vertical cracks 45 . The Brown-Korringa formula applies to the fluid substitution problem in anisotropic rocks and allows to calculate how the compliance of an anisotropic rock changes after saturation with fluid 46 . It assumes low-frequency conditions (pore pressure can equilibrate locally during the passing of seismic wave) which means it is applicable with in-situ seismic frequencies 58 . Hudson's theory describes elastodynamic properties of cracked solids telling how much a distribution of pennyshaped cracks influences the propagation of elastic waves. The wavelengths involved are assumed to be large compared to the size of the cracks and their separation distances.
To track the influence of crack-filling fluid substitution on the SWS we compare the magnitude of S-wave anisotropy of cracked dry ice against cracked watersaturated ice. For this purpose, we use a low-pressure anisotropy background ice model that was inverted for low meltwater time periods (Fig. 4a). This model contains 0.9% of cracks and is the only one for which partial draining of cracks can be considered regarding the meltwater cycle.
Next, we apply the fluid substitution formula 46 to saturate dry cracks with water and obtain stiffness tensor of the water-saturated cracked medium. Knowing the stiffness tensor allows calculation of the magnitude of S-wave anisotropy along various directions (for computation we use the MSAT toolbox 59 ) and comparison against dry case results.
The results are presented in Supplementary Fig. S5. We conclude that for ray paths available in our field data substituting air with water causes a decrease of S-wave anisotropy magnitude. The decrease expected at model space sampled by field data span from −0.15% to 0.1% with a mean of 0.02%, and hence, is an order of magnitude weaker than the influence of cracks itself.
Ice fabric inversion. Methodology for ice fabric description proposed here is a fusion of two methodologies established for describing two different media. The Antarctic ice applications 15,16,44 assume that anisotropic fabric is built from various ice crystals 49 . The shale hydrocarbon reservoirs applications 48 represent the anisotropic fabric by transversely anisotropic shales that may be overprinted with sets of aligned dry cracks 45 resulting in orthorhombic anisotropy. Cracks are introduced using the additional compliance approach 47 .
Here we bring together these two approaches. We assume anisotropic fabric built from various ice crystals 49 and containing vertical water-saturated cracks. Following Verdon et al. 48 we impose dry vertical cracks 45 using the additional compliance approach 47 and then saturate them with water by applying the Brown-Korringa fluid substitution formula 46 . We simulate a low-frequency scenario when there is time for wave-induced crack-pressure increments to flow and equilibrate locally 58 . Similar to Verdon et al. 48 we characterise cracks by crack density and strike only, which enables us to focus on the first order crack density parameter and removes the need to consider other parameters such as crack aperture and rock permeability, to which SWS is less sensitive.
For each icequake cluster-receiver pair median values of the magnitude of S-wave anisotropy and fast S-wave polarization angles are taken. The objective function ζ is based on the L2 norm differences between the observed (obs) and calculated (cal) S-wave anisotropy A and fast S-wave polarization angles ϕ divided by their standard deviations: where n is the number of measurements.

Data availability
Seismometer data from stations of the local seismological network 4D