A Method for Calculating Attenuation in Creeping Materials

Abstract

The phase lag between an applied forcing and a response to that forcing is a fundamental parameter in geophysical signal processing. For solid deforming materials, the phase lag between an oscillatory applied stress and the resulting strain response encapsulates information about the dynamical behavior of materials and attenuation. The phase lag is not directly measured and must be extracted through multiple steps by carefully comparing two time-series signals. The extracted value of the phase lag is highly sensitive to the analysis method, and often there are no comparable values to increase confidence in the calculated results. In this study, we propose a method for extracting the phase lag between two signals when either one or both include an underlying nonlinear trend, which is very common when measuring attenuation in creeping materials. We demonstrate the robustness of the method by analyzing artificial signals with known phases and quantifying their absolute and relative errors. We apply the method to two experimental datasets and compare our results with those of previous studies.

arXiv

Department of Earth and Environmental Science, University of Pennsylvania, PA, USA

Ron Maorronmaor@sas.upenn.edu

We provide a method to accurately extract phase lag between complex signals

The method is evaluated and tested with data from attenuation experiments

We discuss the relevance of this method for studying energy dissipation in geophysical systems and the improved reproducibility in attenuation experiments

Plain Language Summary

Many natural systems are influenced by regular forcing cycles, and their responses to these cycles often show a delay called the phase lag. This delay happens because several mechanisms respond to the cyclic forcing, and some dissipate energy. For geological materials, measuring the phase lag at different forcing frequencies helps to understand the properties of the material and which dynamical processes are involved in the response. This is done by creating an attenuation spectrum. However, since several mechanisms may act at once, it can be challenging to untangle the phase lag caused by a specific forcing frequency. We propose a new method to accurately measure this phase lag, even when multiple processes are involved in the response of materials to forced oscillations. We demonstrate the method’s applicability to both simulated and real laboratory data and provide an open-source algorithm for reproducing our results and using the method in future studies.

1 Introduction

The viscoelastic response of materials to external forcing includes an immediate elastic response and a time-dependent viscous response. As such, the viscous response may be a combination of multiple mechanisms operating simultaneously, each over a different timescale. In geological materials, viscoelastic behavior is tied to numerous phenomena, such as post-glacial isostatic adjustment, tidal dissipation and post-seismic deformation [Karato  Spetzler (1990)]. Deformation in geological materials occurs over long time periods (e.g., 10¹⁵ s for mantle convection), at high homologous temperatures, meaning that these materials exhibit viscous creep in conjunction with viscoelastic behavior [Karato (2008)]. Measuring the viscoelastic response of geological materials in the laboratory is typically done by quantifying the attenuation or internal friction $Q^{-1}$ in forced oscillation experiments performed at elevated pressure and temperature conditions [Jackson  Paterson (1993), Tan . (2001), Cooper (2002), Jackson . (2002), Cao . (2021)]. One of the challenges is that tested samples exhibit both long-term creep and viscoelastic behavior, making data interpretation difficult and susceptible to errors in the estimation of $Q^{-1}$.

Additional complications in estimating $Q^{-1}$ from experiments stem from data analysis techniques that are used to fit experimental results. Attenuation is calculated as the tangent of the phase shift ($\phi$) between two signals in forced oscillation experiments [Cooper (2002), McCarthy  Cooper (2016), Cao . (2021)], typically to the force (or stress) and displacement (strain) signals. Fitting harmonic functions to the stress and strain data entail biases that arise from, for example, assuming that both harmonic signals have the same oscillation frequency, or by detrending a strain signal that also includes a long-term creep component. Even if the spectral separation between the creep behavior and the oscillation frequency is large there is still a problem in fitting multiple parameters for both signals, such as amplitude and phase. Furthermore, because the phase shift of a signal is highly sensitive metric even the length of laboratory-derived data series being fitted to a harmonic function may result in different values of $\phi$.

A common goal of performing attenuation experiments on geological materials is to obtain an attenuation spectrum, which shows the relative energy dissipation in the material at different forcing frequencies. Therefore, there is a need for a robust method to calculate the phase shift at each frequency in a consistent manner, eliminating the additional error inherent in the coefficients of each fit.

We present and test a new method, coined ‘High-Pass Phase Detection’ (HPPD), for calculating the phase shift between two signals (stress and strain) from forced oscillation experiments of creeping materials. This new method employs high-pass signal filtering coupled with the Fast Fourier Transform (FFT) to compute the phase shift through a consistent algorithm, without the need to perform curve-fitting for each signal individually. We show that this new method can consistently recover the phase shift of synthetic signals containing both linear and nonlinear trends. We then apply the new method to experimental data and demonstrate that $Q^{-1}$ can be calculated efficiently and that extracted values are consistent with those of previous studies.

2 Method: High-Pass Phase Detection (HPPD)

Let us consider a general signal composed of a sinusoidal part and an underlying creep trend,

$$x[t]=c(t)+A\sin(2\pi ft+\phi)=\bar{c}[t]+\bar{s}[t],$$ (1)

and its Fourier Transform (FT):

$$X[\omega]=C[\omega]+S[\omega].$$ (2)

The phase of the frequency $\omega$ in the FT spectrum is simultaneously influenced by the sinusoidal component and the creep component. However, if there is enough spectral separation between the primary frequencies constituting the sinusoidal component ($S$) and those constituting the creep component ($C$), filtering methods can be applied to remove the effect of creep on the phase associated with the frequency $\omega$.

To achieve a robust separation between the periodic component that carries the information to estimate $Q^{-1}$ and the creep signal, we use a high-pass response filter with a Butterworth design [Oppenheim  Willsky (1982), Shouran  Elgamli (2020)]. This approach requires specifying a cutoff frequency for the high-pass filter close to the known frequency of the oscillations. Here, we utilize a second-order Butterworth filter. Note that when applying a digital filter, a spurious transient effect occurs in the initial oscillations of the filtered signal [Lynn  Fuerst (1998)]. To address this issue, it is necessary to remove the transient effect from the filtered signal before calculating the phase shift $\phi$.

2.1 Testing with Synthetic Signals

In order to test and quantify the error of the filtering method, we construct synthetic signals with both linear and nonlinear trends and rediscover the phase. The generated signals have the following form:

$$h_{i}(t) =h_{0}\sin(2\pi ft+\phi_{j})+a_{i}t,$$ (3)

$$g_{i}(t) =g_{0}\sin(2\pi ft+\phi_{j})+b_{i}t^{2}+a_{i}t,$$

where $h_{i}(t)$ represents a periodic signal with a linear trend, and $g_{i}(t)$ represents a periodic signal with a nonlinear trend. Both signals have a frequency $f$, a phase $\phi_{j}$, and coefficients $a_{i}$ and $b_{i}$ weighing the trends. We remark that we only deal with the phase and frequency of the signals, hence the units of the other parameters will remain general.

We test the method using phase angles corresponding to the viscoelastic regime ($0<\phi_{j}<\pi/2$), and choose forcing and sampling frequencies that are common in attenuation experiments on geological materials. Figure 1A shows the generated signals that carry the same phase angle. To challenge the method, we deliberately selected a diverse set of trends, spanning a wide range of values that differ by two orders of magnitude.

Figure 1: (A) The synthetic signals from Equation (3). The sampling frequency is 10 Hz, oscillation frequency is 10^-3 Hz, and all the signals have a phase of $\phi=\pi/6$. Dotted lines represent the signals with the nonlinear trend, and solid lines represent the linear trend signals. (B) The same signals after applying the filter. Notice the transient effect at the beginning.

Table 1: Parameters of the signals in Figure 1A and the errors of detecting their phase in Figure 2A. All signals have a frequency of 10^-3 Hz, phase $\phi$  = $\pi/6$, and amplitude $h_{0}=g_{0}=1$. The calibration bias is $\sim$0.756

Signal	$a_{i}$	$b_{i}$	Absolute error [%]
$h_{1}$	$10^{-3}$	-	$1.0\times 10^{-7}$
$h_{2}$	$3\times 10^{-3}$	-	$8.5\times 10^{-5}$
$h_{3}$	$4\times 10^{-3}$	-	$1.27\times 10^{-4}$
$g_{1}$	$10^{-3}$	$10^{-8}$	$8.19\times 10^{-5}$
$g_{2}$	$3\times 10^{-3}$	$10^{-7}$	$7.67\times 10^{-3}$
$g_{3}$	$4\times 10^{-3}$	$10^{-6}$	$4.54\times 10^{-3}$

The underlying trends from all the signals illustrated in Figure 1A are filtered out using a high-pass filter with a second-order Butterworth design. We set the cutoff frequency of the high-pass filter to be half of the oscillation frequency. The results of the filtering are presented in Figure 1B. As shown in the figure, the transient effect of the filter affects the first three oscillations of the signals, which need to be removed from the signals before calculating the phase $\phi$. Because of this spurious transient behavior resulting from the filtering process, having a sufficient number of oscillations in the data is critical when using this method.

Figure 2: (A) The errors between the extracted phase lags for each signal using the filtering method to the known value from constructing the signals. All the errors are below 0.02% error (B) Relative phase lag errors between linear trends, nonlinear trends and a mixture of one linear trend and one nonlinear trend. All the errors are below 0.1% error.

After removing the initial transient oscillations from the filtered signals, we use the Fast Fourier Transform (FFT) to calculate the phase of the dominant frequency in each signal. The generated signals allow us to derive the absolute value of the phase, which is known a priori from constructing the signals in Equation (3) with paramaters summarized in Table LABEL:tab:signals. However, in order to obtain the exact absolute value of the phase, we need to remove the bias introduced by the filter. This bias is a single value attributed to the design of the filter and is consistent across all processed signals. Following this calibration, we calculate the error for each signal, as presented in Figure 2A and summarized in Table LABEL:tab:signals.

When analyzing experimental datasets, the desired parameter is the relative phase angle between the stress and strain signals in each experiment, hence there is no need to calibrate the bias of the filter. To mimic this scenario, we repeat the process described above but now calculate the relative phase angles between different sets of signals: two with linear trends, two with nonlinear trends, and a mixture of a linear trend signal and a nonlinear trend signal. The errors of the relative phase lags are reported in Figure 2B — all of them are below 0.1%.

3 Application: Attenuation Experiments on Halite and PMMA

We test the method using data from attenuation experiments performed with nanoindentation [Badt . (2024)]. This dataset is part of a recent study by [Badt . (2024)] who introduced an experimental technique for measuring the attenuation of halite, olivine, quartz, PMMA and indium by conducting forced oscillation experiments [Cooper (2002), Jackson  Paterson (1993), Cao . (2021)] with a nanoindenter. Figure 3 includes samples of raw data from oscillations at 10^-3 Hz where nonlinear creep trends are present in the strain signals. badt_2024_GRL reported very good agreement between their results and attenuation spectra obtained in previous experiments using different experimental apparatuses. However, in their analysis, badt_2024_GRL estimate the phase shift using a least-squares algorithm, a common approach in previous attenuation studies [Cao . (2021), McCarthy  Cooper (2016), Tan . (2001)].

Figure 3: Raw data from the nanoindentation experiments. The forcing frequency is set to 10^-3 Hz, and the sampling frequency is 1 Hz. $\sigma$ represents the stress with units of milli-Newtons, and $\varepsilon$ represents the strain with units of nanometers. (A) Halite sample. (B) PMMA sample.

Figure 4A shows results from PMMA samples compared to previously published data for PMMA in [Lakes (1998)]. Each blue marker represents a single experiment processed using the method proposed in Section 2. All samples in the figure exhibited a creep trend in the strain signal, while the oscillatory stress signal was smooth and did not require filtering. Each strain signal was filtered using the proposed method, and the transient effect of the filter was removed before calculating the phase between the stress and strain signals. The duration of the transient effect varied among the samples; however, none exceeded four consecutive oscillations. Figure 4B displays the attenuation spectrum calculated for halite. The tests on the halite samples showed a pronounced nonlinear creep trend in all strain signals, but the transient effect did not exceed the length of four oscillations across the spectrum. The deviation of $Q^{-1}$ results in PMMA obtained by nanoindentation from the reported values from [Lakes (1998)] at a frequency $0.1$ Hz (Figure 4A) is due to inherent damping in the nanoindenter [Badt . (2024)].

Figure 4: Log-scale attenuation spectra vs forcing frequency. (A) Attenuation experiments on (A) PMMA and (B) halite using nanoindentation. Blue circles show attenuation results vs frequency calculated using the method introduced in Section 2. Red stars show attenuation results from Lakes (1998) and Badt et al. (2024), respectively.

4 Discussion

Previous studies on attenuation have used nonlinear square fit methods to extract the phase lag [McCarthy  Cooper (2016), Takei . (2014)], often employing the Levenberg–Marquardt algorithm [Levenberg (1944), Pujol (2007)]. This approach requires several steps that are absent in the method presented here. First, there is a need to decide on the functional form of the trend in the signal, such as a linear trend, second-order polynomial, etc. This is usually done by trial and error. Second, the algorithm must be provided with an initial guess for each parameter. This is not straightforward in cases of strong nonlinear creep trends, often occurring at low frequencies, as shown for halite in Figure 4B. Numerical algorithms like Levenberg–Marquardt are very sensitive to initial conditions, adding more complexity and uncertainties to compute the phase lag. In the case of nonlinear trends, there might be six or more parameters to fit, making it difficult to estimate how the error of one parameter is compensated by the others. Moreover, the attenuation is defined as $Q^{-1}=\tan{(\phi)}$, which makes it very sensitive to errors in $\phi$ when values reach $\pi/6$ (i.e., where the linear approximation of $\tan(x)\approx x$ is not valid).

In the method demonstrated here, the signals are processed in a non-iterative way, focusing exclusively on the required phase parameter. This simplifies and strengthen significantly the data analysis process for attenuation experiments and increases the confidence of the calculated values. The only requirement is that the trend of the signal has enough spectral separation from the forced oscillations. This is generally the case since the trend is considered as a very low frequency background signal that does not show a periodic behavior in the entire length of the time series. This fact enables us to detect the phase corresponding to the oscillations of the experiment using FFT accurately. In addition, the process of detecting the phase lag in this study can be easily automated, which is useful when analyzing large datasets or when using numerical models. Rheological models for the upper mantle often involve a power law formulation that leads to nonlinear viscoelasticity [Jain . (2019), Freed  Bürgmann (2004), Masuti . (2023)] and by applying the suggested method, it is possible to use outputs from time-domain power law models to construct detailed attenuation spectra.

Besides attenuation experiments, the phase lag between a periodic forcing and the material response is ubiquitous in fluid systems, for instance, in the circadian radiative forcing of solar radiation in ice-covered waters like lakes. This under-ice radiative heating leads to a periodic yet transient diurnal convective regime over time due to fluid mixing and the nonlinear equation of state of water [Ulloa . (2019)]. The result is a periodic signal in the under-ice water temperature and velocity field with a nonlinear trend, whose actual phase lag response provides information on the transformation rate of available potential energy and kinetic energy in the system [Winters . (2019)]. The proposed method can be readily applied to determine such a phase lag. Other examples include the phase lag between temperature and density in the thermosphere [Kodikara . (2019)] and the oceanic response to wind forcing [Gille . (2001), Zhang . (2021)].

4.1 Conclusions

Currently, there is no consensus on a universally accepted method to determine the ‘phase lag’ $\phi$ between stress and strain signals in attenuation experiments on geological materials (e.g., Cao et. al, 2021). Moreover, the limited amount of data existing at low frequencies ($\leq$ 1 Hz) makes it difficult to assess if the measured $\phi$ and $Q^{-1}$ are consistent. The proposed method (HPPD) for calculating the phase lag presented in this study proves to be effective and relevant for quantifying attenuation in geological materials. We hope the method presented will be useful for future experiments, rendering phase lag extraction more automated and reproducible.

Open Research Section

The data to produce the figures in the article, along with the algorithm to calculate the phase lags are available as Julia code [Bezanson . (2017)] on github.