Rupture Directivity in 3D Inferred From Acoustic Emissions Events in a Mine-Scale Hydraulic Fracturing Experiment

Introduction

Rupture directivity defines preferred earthquake rupture propagation directions through the identification of azimuthal patterns in the apparent seismic source features (e.g., apparent durations, scaled amplitudes), as recorded at a network of stations (Haskell, 1964). Earthquake ruptures can be classified in pure unilateral or pure bilateral (symmetric) ruptures, as well as asymmetric bilateral ruptures, describing an intermediate case between the two previous models. It is well know that a predominance of unilateral ruptures is observed for large earthquakes (McGuire et al., 2002); however, recent studies demonstrated that directivity might be also a common feature of small to moderate events (e.g., Kane et al., 2013; Kurzon et al., 2014; Calderoni et al., 2015; Meng et al., 2020; Ross et al., 2020). So far, rupture directivity has also been identified for a few cases of weak fluid-injection induced earthquakes (Folesky et al., 2016; López-Comino and Cesca, 2018; Király Proag et al., 2019; Wu et al., 2019) and fluid-induced microcracks or acoustic emission (AE) events in mine-scale hydraulic fracturing (HF) experiments (Dahm, 2001), suggesting that such properties could persist at smaller magnitudes. In fact, a recent work has been able to identify directivity effects even at laboratory scale (Kolár et al., 2020).

The most commonly used approach to estimate rupture directivity for small earthquakes is based on empirical Green's functions (EGFs) techniques, using the waveform recordings of fore- or aftershocks with 1–2 magnitude units less than the target event or mainshock (Hartzell, 1978). EGFs are used to account for the poorly know propagation and site effects, which similarly affect EGF and mainshock waveforms, by applying a time-domain or frequency-domain deconvolution process (e.g., Ammon et al., 2006; López-Comino et al., 2012; Abercrombie et al., 2017; Stich et al., 2020). As a result, apparent source time functions (ASTFs) are obtained for each receiver, showing different apparent durations, depending on their relative location to the source and the direction of rupture propagation. Directivity effects and source complexities are typically revealed by the azimuthal patterns of ASTFs in a 2D approach, simplifying the intrinsic 3D rupture propagation problem. Directivity effects in 3D can be also studied applying different approaches, for instance, through stretching techniques (e.g., Warren and Shearer, 2006; Abercrombie et al., 2017; Huang et al., 2017) or using the second order seismic moment (e.g., McGuire, 2004; McGuire, 2017; Fan and McGuire, 2018). The importance of directivity along the dip or vertical direction has been generally evidenced for the analysis of deep-focus earthquakes, by considering the differences among rupture propagation and take off vectors (Warren and Shearer, 2006; Park and Ishii, 2015), as well as for subduction earthquakes (Tilmann et al., 2016; An et al., 2017), using finite source modeling and back projection techniques. For these large scale applications, the directivity analysis is done using surface monitoring networks, i.e. well above the earthquake source, and thus with a poor 3D coverage of the foci. Nevertheless, 3D seismic monitoring configurations deployed recently in deep underground laboratories bear us new opportunities to extend the directivity analysis to microseismicity with moment magnitude (Mw) well below 0 in a 3D environment.

In the last years, decameter-scale in-situ hydraulic stimulation experiments have been carried out in deep underground laboratories with the purpose to improve our understanding about the nucleation and rupture growth processes linked with fluid-injection induced seismicity (Zang et al., 2017; Gischig et al., 2018; Kneafsey et al., 2018; Dresen et al., 2019; Schoenball et al., 2020; Villiger et al., 2020; Villiger et al., 2021). These experiments have been promoted in recent years by the development of enhanced geothermal systems and unconventional resources driven by HF operations (Schultz et al., 2020). Mine-scale in-situ experiments serve as a bridge between laboratory experiments and deep reservoir stimulations, offering improved control and more realistic boundary conditions. AE piezoelectric sensors are deployed at different borehole locations around the fluid-injection intervals providing complex 3D seismic monitoring networks. Microfractures generated during HF stimulations are then characterized by AE signals involving, for instance, M_w between −4.2 and −3.5 and rupture sizes of decimeter scale (Kwiatek et al., 2018).

In this framework, on June 2015 a mine-scale underground HF experiment at 410 m depth was carried out in the Äspö Hard Rock Laboratory (HRL), Sweden (Zang et al., 2017). Six HF stimulations with three different injection schemes were tested to assess the fatigue HF (FHF) concept in order to reduce the number and magnitude of fluid-injection induced seismicity. Each HF stimulation composed by the initial fracture phase (F) and up to five refracturing phases (RF) was performed at selected injection intervals free of pre-existing fractures along a subhorizontal 28 m-long borehole (Figure 1). This injection borehole was drilled in the orientation of the minimum horizontal stress magnitude (S_h) with a strike of ∼N210 W. The 3D AE monitoring network consists of 11 AE uniaxial side view sensors, recording in the frequency range of 1–100 kHz, oriented toward the stimulated volume of 30 × 30 × 30 m and deployed in different boreholes and along existing experimental tunnels of Äspö HRL. These AE sensors are entirely based on the piezoelectric effect, which means that the seismic wave is guided directly into the sensor where it generates an electric output signal proportional to stress changes introduced. The data acquisition system was improved to operate with 1 MHz sampling rate, obtaining high-quality full waveforms. The AE activity generated during each HF stage has been well-studied and characterized using both triggered mode (Zang et al., 2017; Kwiatek et al., 2018) and continuous recordings (López-Comino et al., 2017; Niemz et al., 2020). Figure 1 shows the largest AE events with magnitudes (relative magnitude, M_AE) exceeding 2.8, according to the last updated catalog, which contains 4,302 events with M_AE between ∼2.0 and ∼4.0 (Niemz et al., 2020). This enhanced catalogue incorporating weaker events allowed for the identification of planar seismogenic regions associated to each HF stimulation (Niemz et al., 2020). On the other hand, the M_w range has been estimated between −4.2 and −3.5 for the 196 AE events recorded in the triggered mode; furthermore, 21 moment tensor inversions were performed, identifying heterogeneous focal mechanisms and suggesting a complex microfracture network generated for the hydraulic stimulations (Kwiatek et al., 2018). Interestingly, for the large majority of the retrieved focal mechanisms, none of the potential fault planes matches the orientation of the larger scale planar seismogenic regions.

FIGURE 1. (A) and (B) AE hypocenters with M_AE >2.8 located in the hydraulic fracturing experiments in the Äspö HRL from the last updated catalogue of Niemz et al. (2020). Stars denote the largest AE events analyzed in this study, indicating the corresponding M_AE. The AE sensors (triangles) and injection intervals (squares) along the 28-m-long injection borehole (blue line) are also shown. The direction of view in (B) is indicated in (A). (C) Temporal evolution of the AE events of each HF experiment shown in (A) is plotted together with the injection pressures (black lines) and flow rates (blue lines). Initial time (t₀) is indicated for each plot. HF1, HF2 and HF6 were performed with a conventional continuous injection scheme, while a cyclic progressive injection scheme is followed in HF3 showing less AE activity and lower M_AE.

Typically, the source finiteness of small magnitude earthquakes, as well as AE events, is neglected under a point source approximation. Seismic source properties are then reduced to earthquake location, magnitude, and, in few cases, to a moment tensor, which analysis has been based on P wave first-motion polarities (Kwiatek et al., 2016). Relevant rupture features, such as rupture size, duration and directivity are not discussed and still represent a challenge at the scale of microseismicity and AEs. Here, we characterize the source finiteness for the largest AE events recorded during the HF experiments in the Äspö HRL, through a novel 3D-directivity analysis applied to 3D seismic monitoring networks. Unilateral ruptures, as well as, symmetric and asymmetric bilateral ruptures are tested to find the best fitting model for the 3D pattern of ASTFs. Such advanced characterization of the rupture processes can be used to discriminate the rupture and auxiliary plane and to better understand the geometry of induced fractures.

Methods

Apparent Durations in Mine-Scale Experiments

The high resolution AE seismicity recorded during the different HF stimulations in combination with the 3D geometry of the monitoring network and the high sampling of recordings at the Äspö HRL provide us a great opportunity to obtain apparent durations for AE events, which is a challenging step, considering their extremely low magnitudes and short durations. For some cases of global seismological applications, for instance in the case of deep-focus earthquakes, the initial portion of P-wave displacement is isolated from most other phases and represents well the source time function with a scalar correction for the seismic moment (Fukao, 1972; Kikuchi and lshida, 1993; Beck et al., 1995; Tibi et al., 1999). Similar conditions are found for the first arrivals of the P-waves recorded in the AE network of this small-scale experiment. The average distance from the source to the borehole sensors is only about 18 m, and the wave propagation is well described by a homogeneous full-space (Niemz et al., 2020). Therefore, the effects of local structure are minimal and secondary arrivals, for instance from reflections at the gallery walls, are expected to be either weak or not significantly overlapping with the fairly short duration of the direct P wave. In such circumstances, we can assume the duration of the first pulse observed in the P-wave as a rough estimation of the apparent duration. Although the amplitudes and polarities of these waveforms are affected by the source radiation pattern of the moment tensor and eventually by the sensor orientation and coupling (Manthei et al., 2001; Maghsoudi et al., 2013), such effects do not alter the apparent durations, which are only considered for the rupture directivity modeling.

We manually select six AE events, four from HF2 and two from HF6 (Figure 1), where the first P-wave pulses of unfiltered waveforms are well identified in at least nine out of 11 AE sensors (Figure 2, Supplementary Figure S1). Normalized square amplitudes of the first P-wave pulses are plotted to improve the visualization of the apparent duration measurements. The M_AE of these selected AE events range from 3.49 to 3.92 (i.e., M_w ∼ −3.5, according to Kwiatek et al., 2018) corresponding with the largest and high-quality AE events. Note that few AE signals are discarded from our analysis, because the first P-wave pulse is not well-constrained, in presence of noisy records, low P amplitudes (e.g., for nodal stations affected by the radiation pattern near the intersection between the compressional and dilation quadrants) or anomalous long pulses, compared to the remaining stations (Figure 2). P-wave arrivals are manually picked and we calculate automatically the end of the apparent duration at the intersection of the first P-wave pulse with the baseline. Some variations in the apparent durations are observed in our target AE events ranging from 0.016 to 0.053 ms. This implies apparent corner frequencies in the range of 18–62 kHz and magnitudes as small as M_W ∼ −3.5 (Kwiatek et al., 2018). On the other hand, the frequency dependent sensitivity of the AE sensors is known to affect recordings especially above 100 kHz (Ono 2018) and include low frequency noise below 3 kHz (Niemz et al., 2020), which could potentially affect our estimate of the apparent durations. To confirm the robustness of our results, we compared our apparent duration estimates, as retrieved from data filters in the frequency bands 0.5–500 kHz, 0.5–100 kHz, and 3–50 kHz (Supplementary Figure S2), and only find out negligible differences yielding the same apparent duration measurements. The lowest possible upper bound of the previous applied bandpass filter can be defined in 31 kHz because the first P-wave pulse cannot contain lower frequencies than the inverse of its period (about two times the minimum apparent duration). The resolution of finite source patterns, such as the directivity, formally requires the inclusion of higher frequencies (e.g., above 50 kHz). While apparent duration results from 2D directivity analysis are typically shown sorting based on the station azimuth, here we introduce other sorting strategy to visualize the 3D directivity effects based on the parameter κ described in the next section.

FIGURE 2. Apparent durations for the largest AE event recorded during HF2 in the Äspö HRL. We plot the unfiltered waveforms (black lines), P-wave window (orange band), the first P-wave pulses (red lines) and the apparent durations (dots according to the colorbar in Figure 3A). Manual picking and the baseline are indicated in the second column with vertical dashed gray lines and horizontal gray lines respectively. Waveforms are normalized according to the maximum amplitude (see labels in each trace of the first column) and sorted by the parameter κ (see labels of the second column). Some AE signals are discarded from the rupture directivity analysis because the first P-wave pulse is not well constrained and does not allow a clear estimate of the apparent durations (grey traces in second column for AE07 and AE08). The third column shows the normalized square amplitude of the P-wave pulse shown in the second column to improve the visualization of the apparent duration measurements.

Rupture Directivity Analysis in 3D

The source parameters associated to a line source (Haskell, 1964) are commonly inferred through the apparent duration fitting (e.g., Cesca et al., 2011; López-Comino et al., 2016). For unilateral ruptures, the apparent durations, Δt(ϕ), show the longest duration in the forward direction of the rupture propagation (α) and the shortest duration in the backward direction, which can be written as:

depending on the azimuth (ϕ) between the source and each receiver, and involving different variables: rise time (t_r), rupture length (L), rupture velocity (v_r), and P (or S) wave velocities (v_P,S). In this manner, the rupture directivity is constrained along a horizontal plane and defined by a single angle, α.

For our 3D mine-scale case, we will consider spherical coordinates involving two angles (ϕ, θ): the azimuth ϕ, previously described in Eq. 1, and a polar angle θ, describing the deviation from the vertical axis (Supplementary Figure S3). We define the 3D rupture directivity vector, $\overrightarrow{d}\left(\varphi ,\theta \right)$, applied at the source location, defined by these two angles, and, in the same way, the station vector, $\overrightarrow{s}\left({\varphi }_i,{\theta }_i\right)$, from the source location to each station (i). Thus, we introduce the parameter κ by:

representing the angle between the directivity vector and the station vector, ranging from 0° to 180° (Supplementary Figure S3). For κ values near 0°, indicating an alignment among station and rupture vectors, and thus observations ahead of the rupture direction, we expect to observe the shortest apparent durations. Conversely, we expect the opposite behavior, showing the largest apparent durations, for κ-values of ∼180°.

Replacing the angle ϕ-α in Eq. 1 with the new parameter κ, we can model the apparent durations Δt in 3D as a function of two angles, i.e., Δt (ϕ, θ), assuming a unilateral rupture. Consequently, the general 3D case for unilateral and (asymmetric) bilateral ruptures can be defined combining κ and the parameterization of Cesca et al. (2011), by the equation:

where the asymmetry of the rupture (χ) ranges from 0 for pure unilateral rupture to 0.5 for pure (symmetric) bilateral rupture. Theoretical radiation patterns of apparent durations for three significant cases (pure unilateral, pure bilateral and asymmetric bilateral) are shown in Supplementary Figure S4 for one specific directivity vector and considering reasonable values for the involved parameters according to the results obtained in the next section. Average apparent durations increase with decreasing rupture velocities, increasing rupture length and increasing rise time.

Theoretical curves of apparent durations vs. κ values following Eq. 3 can then be tested for all directions, applying a full grid search for ϕ in [0°, 360°] and θ in [0°, 180°]. L1-misfit for each tested direction is then obtained by comparing theoretical and observed apparent durations at each station, showing the uncertainties associated to the best solution (minimum L1-misfit) of the directivity vector. The inversion is solved through a full grid search approach including reasonable intervals and increments for the involved parameters. The χ-parameter varies between 0 (pure unilateral rupture) to 0.5 (pure bilateral rupture). Average P-wave and S-wave velocities were measured in-situ (v_p = 5,800 m/s and v_s = 3,200 m/s, Zang et al., 2017), thus we can consider a rupture velocity interval between 0.7 v_s and 1.0 v_s. The rise time is not expected to be larger than 1/3 of the true duration (e.g., Stein and Wysession, 2003), thereby we test below 1/3 of the average apparent duration. Despite our knowledge about the rupture length is limited for this magnitude range, some reasonable values of 10–30 cm can be estimated for such range of magnitudes (e.g., Eshelby, 1957; Kwiatek et al., 2018); then, a broader interval from 1 to 30 cm is tested. A jackknife test is also applied to assess the uncertainties of each parameter, which are quantified by the mean and standard deviation calculated from the different iterations (see Supplementary Table S1).

Directivity of Decimeter-Scale Acoustic Emission Events

Theoretical predictions assuming unilateral and asymmetric bilateral rupture models are considered to adjust the 3D radiation pattern defined by the apparent durations previously identified at each AE sensor (Figure 3A). Such predictions are also compared using the 2D approach yielding similar results and azimuthal directivity directions (Supplementary Figure S5). We scan all possible directions in 3D by increments of 5° for both vector angles (ϕ, θ) and the L1-misfit is calculated for each tested direction (Figure 3B, Supplementary Movies S1 and S2). The resulting rupture directivity vector is then defined by the minimum L1-misfit. Figures 3, 4 summarize the resulting inverted parameters from our 3D rupture directivity analysis (see also Supplementary Table S1 and Supplementary Movie S3).

FIGURE 3. Rupture directivity analysis in 3D for the target AE events recorded during HF2 and HF6 in the Äspö HRL. (A) For each target AE event, we show the apparent durations measurements (dots) along with the theoretical prediction (brown line) for the inverted model using a full grid search approach (labels indicate the solution of each parameter). (B) L1-misfit for each tested direction (dots) for the largest AE events recorded during HF2 (top) and HF6 (bottom). The directivity vector solution with minimum L1-misfit (red arrow) and the station vectors (gray lines) for each AE sensor are also shown. Red dots indicate the origin of the coordinate system. The radius for each tested direction is also scaled according to the L1-misfit values as in the colorbar.

FIGURE 4. 3D views of the resulting rupture directivity vectors for the four largest target AE events in the Äspö HRL: (A) top view and (B) perspective view. Planar seismogenic regions identified by the AE activity (Niemz et al., 2020) are shown considering the hypocentral locations of our target events. The AE sensors are located in nearby boreholes and tunnel; they are oriented toward the stimulated volume along the 28-m-long injection borehole. Focal mechanisms obtained by Kwiatek et al. (2018) for the largest AE event in HF2 and HF6 are shown with their rupture directivity vectors (brown arrows) in (A). M_AE is indicated for the target events in (B).

Despite our monitoring network configuration includes only AE sensors located above the injection borehole, robust results are inferred for four out of six target AE events (M_AE > 3.52), showing the best adjustments for L1-misfits < 0.003 ms and small uncertainties of <10° associated to the 3D-directivity directions defined by the ϕ and θ angles. These results reveal a slightly predominance of unilateral ruptures for three out of four AE events. Similar azimuthal directivity directions are observed for the AE events belonging to the same HF stimulations, being 200°–255° for HF2 and 20°–40° for HF6 (Figures 3A, 4A). A predominant pattern along the vertical direction cannot be determined. Beside some intrinsic trade-offs among the involved parameters (Stich et al., 2020; López-Comino et al., 2021), we clearly identify high rupture velocities (3.14–3.20 km/s), very short rise times (<0.002 ms), and rupture lengths ranging 10–18 cm.

Discussion and Conclusions

Rupture directivity analysis is extended at mine-scale and can be recovered for microearthquakes recorded by 3D complex seismic monitoring configurations. The high-quality full waveforms recorded for the largest AE events during the HF2 and HF6 stimulations in the Äspö HRL has allowed us to estimate the apparent durations through the first pulses observed in the P-wave of the AE signals. We described a 3D-directivity approach where theoretical radiation patterns of apparent durations for unilateral and (asymmetric) bilateral ruptures are characterized by the introduction of a parameter κ, representing the angle between the directivity vector and the station vector. A full grid search is proposed to solve the inversion problem assessing the model performance through the L1-misfit in all plausible rupture directions and quantifying the associated uncertainties by a jackknife approach.

We recognize, for the first time, directivity effects for microfractures within granitic rock with magnitude as small as M_w ∼ −3.5 and rupture length of ∼13 cm lengths. Given the small number and narrow band of target magnitudes, no empirical relationships between magnitude and rupture length can be derived from our results (Supplementary Figure S6). Similar average apparent durations $\left(\overline{t}\right)$ are observed for all target AE events. Thereby, in these cases the rupture length is mostly controlled by the rupture propagation mode, being larger for (asymmetric) bilateral ruptures and shorter for unilateral ruptures. Accordingly, when considering fixed values of rupture length, rise time and rupture velocity, the apparent durations decrease from unilateral to pure bilateral ruptures (Supplementary Figure S4).

Our results can resolve the fault plane ambiguities for the largest AE events in HF2 and HF6, when taking into account the moment tensor solutions reported in Kwiatek et al. (2018). Indeed, rupture vectors are always coplanar to one of the potential fault planes (see Figure 4A and the indicated rupture directivity vectors in the moment tensor solutions). Similar conclusions were drawn from a moment tensor study of hydrofrac-induced AE events in salt (Dahm et al., 1999). On the other hand, multiple fractures planes were inferry a 2D approximation of the ellipsoidal AE clusters revealing some strike variations and steep dipping angles (Niemz et al., 2020). Rupture directivity vectors can be discussed with respect to these main planes, considering the hypocentral locations of target AE events (Figure 4). For example, a secondary fracture zone mapped in HF2–RF4 and HF2–RF5, with a predominantly grew during the shut-in phase, is roughly aligned with the rupture vector of the largest AE event, which is controlled by the interaction with the previously fractured zone of HF1. For the remaining studied AEs, however, the rupture planes that we identify are inclined, when not perpendicular, to the larger scale (few meters) planar seismogenic regions(Figure 4B). Co-seismic rupture planes out of the plane of the hydrofracture were also found by a directivity analysis of AE in salt rock (Dahm, 2001). This confirms a clear difference among the orientation of the large scale extension of the fractured region, and the small scale orientation of single fracture planes. Indeed, single focal mechanisms already indicated a significant variability (Kwiatek et al., 2018), except for their pressure axes being consistently oriented according to a common stress field (i.e., σ_H oriented NE-SW to NNE-SSW).

The discussed mine-scale HF experiments in deep underground laboratories bear us also the chance to discuss rupture directivity patterns, which have been proposed for larger scale seismicity. Fluid-induced seismic processes have been observed over a broad range of spatial scales (Davies et al., 2013; Schultz et al., 2020), and very often with dense dedicated networks, which allow to resolve the rupture propagation for small to moderate earthquakes (e.g., 1.0 < M_w < 5.8). Lui and Huang (2019) calculated the directivity effects of the largest induced Oklahoma earthquakes and attributed the difference in rupture directions to expected pressurization of the fault zone, which relates to the distance away from injection zones and total injected volume. Folesky et al. (2016) analyzed the rupture directivity of the largest seismic events associated with the stimulation of geothermal reservoir in Basel (Switzerland) and found that the preferred rupture propagation depends on magnitude; events with local magnitude larger than two propagated backward into the perturbed volume while smaller events propagated away from the well. These studies suggest that the fluid injection has an effect on directional properties of the earthquake rupture processes. In this sense, our analysis evidences a predominant pattern of rupture directivity directions, as they are in general well aligned with the orientation of the minimum horizontal stress (e.g., Nolen-Hoeksema and Ruff 2001; Zang and Stephansson 2010); a preferred rupture directivity toward or away from the injection point is, however, not identified.

While based on a limited dataset, our findings suggest that the local stress condition may control the predominant direction of the rupture growth, beside the large scale distribution of AE locations and the pressure axis of the focal mechanism. This study confirms that directivity effects persist even at decimeter-scale ruptures and extends our knowledge to better understand the triggering processes of fluid-injection-induced earthquakes.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The datasets selected for this study are available via request to the authors. Catalogue of induced AEs from Niemz et al. (2020). Moment tensor solutions were extracted from Kwiatek et al. (2018).

Author Contributions

JL-C and SC planned the research. JL-C has realized the rupture directivity analysis, supervised by SC; PN compiled the data and realized the movies of the results. JL-C prepared the initial draft and figures, with contributions and editing from all authors. AZ was the principal investigator of the 2015 mine-scale Äspö HRL underground experiments. All authors have contributed to the interpretation and discussion.

Funding

This research has been funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N° 754446 and UGR Research and Knowledge Transfer Found – Athenea3i; and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer (407141557). SC and PN received funding from the EU project PostMinQuake, RFCS (899192). AZ was supported by funding received from the European Union’s Horizon 2020 research and innovation programme, Grant Agreement No. 691728 (Destress).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor (APR) declared a past co-authorship with one of the authors (TD).

Acknowledgments

The in situ experiment at Äspö Hard Rock Laboratory (HRL) was supported by the GFZ German Research Center for Geosciences (75%), the KIT Karlsruhe Institute of Technology (15%) and the Nova Center for University Studies, Research and Development Oskarshamn (10%). An additional in-kind contribution of the Swedish Nuclear Fuel and Waste Management Co. (SKB) for using Äspö HRL as test site for geothermal research is greatly acknowledged.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.670757/full#supplementary-material

References