Effects of microstructural heterogeneity on the macroscopic spectrum of elastically accommodated grain-boundary sliding

To validate the numerical framework, we first computed the macroscopic response of a perfectly regular hexagonal tessellation (Figure 2, dashed black lines). Here and in all later simulations, a Poisson’s ratio of the grain interiors $\nu=0.35$ was adopted as in \citeAghahremaniEffectGrainBoundary1980. The resulting spectra reproduce the analytical and numerical results of \citeAghahremaniEffectGrainBoundary1980. In this idealized geometry, all grain boundaries have identical lengths and all triple junctions meet at symmetric $120^{\circ}$ angles, strongly constraining intergranular sliding. As a result, the aggregate exhibits a relatively high relaxed storage modulus ($G^{\prime}/\mu\approx 0.82$ as $\omega\rightarrow 0$) and a single narrow Debye-like attenuation peak with $G^{\prime\prime}_{\max}/\mu\approx 0.09$. It is common to define the relaxation strength $\Delta$ as:

In this case the calculated relaxation strength is $\Delta=0.18$.

We then replaced this idealized benchmark with irregular polygonal grains generated from a Voronoi tessellation with a narrow grain-size distribution ($\sigma_{a}=0.05$). As shown in Figure 2a, the introduction of randomized boundary lengths and irregular junction angles allows greater intergranular sliding, leading to a larger relaxed compliance. Consequently, the relaxation strength rises sharply to $\Delta\approx 0.28$.

In addition to increasing the total relaxation strength, the transition from regular hexagons to irregular polygonal grains systematically shifts the attenuation peak to lower frequency by a factor of approximately $1.4$. This offset remains largely unchanged as $\sigma_{a}$ is increased in the subsequent simulations, indicating that it is not controlled by the width of the grain-size distribution itself. It therefore appears to arise mainly from grain-shape and network-topology irregularity relative to the regular hexagonal benchmark.

Having established the strong effect of irregular grain geometry, we next investigate the role of grain-size variance within the Voronoi aggregates by increasing the standard deviation of the log-normal grain-size distribution from $\sigma_{a}=0.05$ to $\sigma_{a}=0.50$.

As shown in Figure 2, broadening the grain-size distribution produces only a modest change in the macroscopic spectra. In Figure 2a, increasing $\sigma_{a}$ produces a slight increase in the relaxed modulus. Correspondingly, Figure 2b shows a minor but systematic reduction in the peak attenuation height. However, the loss spectrum remains narrow and retains a Debye-like form across the full range of $\sigma_{a}$ examined.

Increasing grain-size variance does not produce appreciable spectral broadening, nor does it significantly alter the frequency shift introduced by the transition from regular to irregular grain geometry. Thus, although grain-size heterogeneity slightly modifies the relaxation strength and peak position, it does not broaden or suppress the attenuation peak to account for the weak or poorly resolved peak observed in dry olivine experiments. Within the present model framework, geometric heterogeneity alone is therefore insufficient to explain the experimental discrepancy.

To test this hypothesis, we assigned each internal grain boundary an independently sampled viscosity factor in our 2-D Voronoi aggregates, drawn from a log-normal distribution: $\ln\eta_{\mathrm{gb}}\sim\mathcal{N}(\bar{\eta},\sigma_{\eta}^{2})$. The geometric mean viscosity $\bar{\eta}$ was held fixed so that the characteristic relaxation frequency remained within the numerical frequency window, and the distribution width ($\sigma_{\eta}$) was varied from 0.10 (near-uniform) to 7.00 (highly heterogeneous).

The resulting macroscopic attenuation spectra, computed as the ensemble-averaged loss modulus $G^{\prime\prime}/\mu$, are presented in Figure 3a. As the variance in boundary viscosity increases, the localized Debye-like peak systematically flattens, progressively broadening into a weaker, more distributed background. Figure 3b quantifies this transition: as $\sigma_{\eta}$ widens, the peak height decreases monotonically while the full-width at half-maximum (FWHM) increases approximately linearly.

To examine the origin of this spectral broadening, we grouped grain boundaries into ten log-spaced bins according to their viscosities and decomposed the total macroscopic dissipation into additive contributions from each bin (Figure 4). Details of the decomposition procedure are given in D.

At near-uniform distributions ($\sigma_{\eta}=0.10$), all viscosity bins relax nearly simultaneously, superimposing to form a single, narrow Debye-like peak. However, as $\sigma_{\eta}$ increases, the individual bin responses separate along the frequency axis. Each viscosity class retains a localized Debye-like response, and the broad macroscopic spectrum arises from their superposition over distinct characteristic frequencies.

Although the full 2-D simulations resolve a mechanically coupled network of grains, boundaries, and triple junctions, it is useful to determine whether the integrated dissipation is controlled primarily by network connectivity or by the total abundance of each viscosity class.

To assess this, we integrated each bin-resolved loss spectrum over frequency and compared the resulting dissipation with the total boundary length in that viscosity class. We then compared this integrated dissipation against the corresponding total boundary length of that viscosity class within the aggregate (Figure 5). For most distribution widths, the two quantities exhibit an approximately linear relationship.

This linearity implies that, to leading order, the macroscopic dissipation contributed by a given viscosity class scales with its total fractional representation in the grain-boundary network. In other words, although the instantaneous stress and sliding fields remain fully coupled through the elastic grain framework, the integrated dissipation can be partitioned approximately additively across viscosity classes. The dominant control on the total energy dissipated by a given class is therefore its total boundary length, rather than higher-order details of boundary connectivity.

Small deviations from linearity become more apparent only at the largest distribution widths ($\sigma_{\eta}\gtrsim 5$). In these cases, the characteristic relaxation frequencies of the most viscous and least viscous bins are shifted toward or beyond the edges of the numerical frequency window, so their integrated areas are slightly underestimated. These departures are therefore attributable to finite frequency-window truncation rather than a signature of non-linearity.

This result has two important implications. First, it shows that broad aggregate-scale attenuation can be interpreted as the cumulative contribution of many localized relaxations associated with different grain-boundary viscosities. Second, it provides the numerical motivation for a reduced-order rheological model in which the macroscopic EAGBS spectrum is approximated by a weighted superposition of localized Debye-like responses.

Figure 5 shows that the integrated dissipation associated with each grain-boundary viscosity class scales approximately with its fractional boundary length in the aggregate. This suggests that, to leading order, the macroscopic EAGBS response can be approximated as a sum of contributions from different viscosity classes.

This motivates a reduced-order (0-D) representation of the EAGBS spectrum. In this approximation, the macroscopic response is represented as the superposition of many Debye relaxation elements (standard linear solids) connected in parallel, each associated with a grain-boundary viscosity $\eta_{\mathrm{gb}}$. In continuous form, the loss modulus is

where $g^{\prime\prime}_{\mathrm{single}}(\omega;\eta_{\mathrm{gb}})$ is the localized Debye-like loss response of an individual viscosity class and $p(\eta_{\mathrm{gb}})$ is the corresponding length-weighted probability density function. In practice, the reduced-order model is implemented at the level of the complex modulus, so that both the storage and loss components, $G^{\prime}(\omega)$ and $G^{\prime\prime}(\omega)$, are obtained from the real and imaginary parts of the same distributed standard-linear-solid representation. The explicit form used for this calculation is given in Appendix E.3.

Figure 6 shows that this simple model reproduces the main qualitative trend seen in the FEM calculations: increasing the width of the viscosity distribution suppresses the peak height and broadens the spectrum into a weak background.

The 0-D formulation provides both a compact physical interpretation of the numerical results and a computationally efficient way to explore how an assumed grain-boundary viscosity distribution maps into an attenuation spectrum.

The 0-D model does not eliminate the role of geometry. As shown in Section 3, irregular grain topology changes both the relaxation strength and the characteristic frequency relative to the regular hexagonal benchmark. In the narrow-distribution limit, these geometric effects provide the calibration for the amplitude and frequency scale of the reduced-order model. The present formulation should therefore be understood as separating the problem into two parts: geometric calibration of the relaxation response, and statistical averaging over the grain-boundary viscosity distribution.

Figure 7 shows that dry olivine is characterized by a broad attenuation background with only a weak peak-like feature expressed mainly as a change in slope [quOnsetAnelasticBehavior2024]. The results of Sections 4 and 5.1 suggest a possible explanation for this behaviour. If grain boundaries sample a sufficiently broad distribution of viscosities, then the EAGBS contribution is spread over a wide range of relaxation times and may become difficult to distinguish from the background attenuation.

Such a broad distribution is physically plausible in dry olivine. Grain-boundary character in olivine spans a large five-dimensional crystallographic space [marquardtStructureCompositionOlivine2018]. Since olivine is orthorhombic with low crystal symmetry, its grain-boundary character space is far more diverse than that of face-centered cubic crystals such as aluminum. SEM observations further suggest that olivine grain boundaries do not cluster strongly around a small set of preferred misorientation axes [faulGrainMisorientationsPartially1999]. Molecular-dynamics calculations indicate that grain-boundary transport properties vary by more than three orders of magnitude across this space [wagnerAnisotropySelfdiffusionForsterite2016]. Additional variability in boundary thickness, local defect content, and interfacial structure would broaden the effective viscosity distribution still further. Dry olivine is therefore expected to sample a wide range of grain-boundary viscosities.

For scale, a log-normal distribution with $\sigma_{\eta}=4$ has a one-standard-deviation interval of $\exp(\pm 4)$ about the median, spanning approximately $3\times 10^{3}$, or about 3.5 orders of magnitude, in viscosity. This is comparable to the order-of-magnitude variability suggested by molecular-dynamics calculations, although the precise width of the effective grain-boundary viscosity distribution in dry olivine remains poorly constrained.

The same framework may also explain why pronounced EAGBS peaks are more readily observed in simpler polycrystalline systems. In face-centered cubic metals such as aluminum, systematic surveys of all 388 distinct grain-boundary types find that the spread in mobility is narrower than in olivine [olmstedSurveyComputedGrain2009]. The sliding time is then concentrated around a more localized characteristic timescale, allowing the EAGBS peak to survive macroscopic averaging. This is consistent with classic observations in metals [keExperimentalEvidenceViscous1947]. The contrast with dry olivine is therefore consistent with the broader conclusion of this study: whether EAGBS appears experimentally as a distinct peak or only as a weak background contribution depends strongly on the width and position of the underlying grain-boundary viscosity distribution.

Wet olivine provides an instructive contrast. In Figure 7, the wet spectra show a more localized and higher-amplitude peak-like feature than the dry spectra. Liu et al. (2023) interpreted this feature in terms of bulk diffusion of hydrogen defects rather than grain-boundary sliding, mainly because the fitted peak position showed little clear dependence on grain size in their analysis [liuEffectWaterSeismic2023]. Karato (2025), however, argued that this conclusion is not robust and that, after excluding an outlier, the grain-size dependence is more consistent with an EAGBS interpretation [karatoCausalityItsImplications2025]. The present results suggest one possible reason why an EAGBS peak may be more visible in wet olivine: hydration may modify the grain-boundary state and thereby narrow the effective viscosity distribution.

Atomistic calculations on forsterite show that hydrogen-bearing defects are energetically favored at grain boundaries and that hydrated boundaries may develop lower-density interfacial regions [deleeuwProtoncontainingDefectsForsterite2000]. This suggests that hydration may modify grain boundaries before comparable effects are expressed uniformly throughout the bulk. If hydration drives boundaries toward a more liquid-like or weakly structured interfacial state, then the broad dry-olivine distribution of $\eta_{\mathrm{gb}}$ may narrow because the effective film viscosity would depend less strongly on the original crystallographic boundary character. A qualitatively similar effect is seen in melt-bearing olivine aggregates, where a more distinct peak-like relaxation feature is observed [faulShearWaveAttenuation2004, jacksonShearWaveAttenuation2004]. Within the present framework, such observations are consistent with the idea that boundary-state modification may make EAGBS easier to observe experimentally.

This interpretation does not exclude Liu et al.’s bulk hydrogen-defect diffusion explanation, but it suggests a simple experimental discriminator. If the wet peak is controlled primarily by a bulk hydrogen-defect relaxation, then its strength should continue to increase with water content. By contrast, if the peak reflects grain-boundary sliding made more visible by preferential boundary hydration and boundary-state modification, then its amplitude should eventually saturate once most grain boundaries have been modified. In that case, further water uptake would produce only limited additional enhancement because the boundary network would already be close to saturation.

Upper-mantle viscoelasticity is expressed through both seismic attenuation and velocity dispersion. For weak attenuation, the inverse quality factor is approximately

and the corresponding shear-wave velocity is

where $\rho$ is density.

Recent work has argued that transient diffusion creep alone may be insufficient to explain the observed attenuation of the asthenosphere at seismic periods. In particular, \citeArudgeViscoelasticRheologyTransient2026 showed that extrapolation of transient diffusion creep to seismic frequencies yields attenuation significantly smaller than the $Q^{-1}\sim 10^{-2}$ often inferred from seismology. This suggests that an additional dissipative mechanism may be required under upper-mantle conditions.

In terms of dissipation, \citeApriestleyThermalAnisotropicStructure2024 parameterized surface-wave velocity as a function of temperature and pressure by combining a low-temperature elastic dependence of density and unrelaxed shear modulus with a high-temperature Andrade-type anelastic reduction. They emphasized that this high-temperature parameterization is empirical, although it provides a simple description of the observed broad power-law anelastic behavior. Here we retain the same low-temperature elastic parameterization, but replace the empirical high-temperature Andrade-type reduction with the distributed-viscosity EAGBS model developed in Section 5.1, to examine whether a physically motivated grain-boundary sliding mechanism can generate a comparably broad anelastic velocity reduction. The reference sliding timescale is anchored to the dry olivine experiments of \citeAquOnsetAnelasticBehavior2024, while the activation energy of grain-boundary relaxation, $E_{a}$, and the width of the log-normal grain-boundary viscosity distribution, $\sigma_{\eta}$, are treated as fitting parameters. A full mathematical statement of the fitted reduced-order model, including the definitions of $\mu(T,P)$, $\rho(T,P)$, and $V_{s,u}(T,P)$, together with the inversion procedure and best-fitting parameters, is provided in  E and Table 1.

The resulting fit and the corresponding attenuation calculated at a 40 s period are shown in Figure 8. The fit shows that the distributed-viscosity EAGBS model can account for the nonlinear high-temperature reduction in $V_{s}$, while also generating attenuation of the order required by asthenospheric observations. In particular, the narrow-distribution case produces insufficient attenuation, whereas the fitted broader distribution, with $\sigma_{\eta}=3.92$, yields attenuation of order $Q^{-1}\sim 10^{-2}$ over asthenospheric temperatures.

The best-fitting activation energy is $E_{a}=429$ kJ/mol, which is slightly lower than, but still comparable to, the range of grain-boundary diffusion activation energies inferred experimentally by \citeAyabeGrainBoundaryDiffusionCreep2020 and 666 kJ/mol by \citeAquOnsetAnelasticBehavior2024. The fitted viscosity-distribution width is $\sigma_{\eta}=3.92$, indicating that substantial grain-boundary rheological heterogeneity is required if EAGBS is to account simultaneously for the observed velocity reduction and seismically relevant attenuation. This fitted value is also broadly consistent with the hypothesis developed in Sections 4 and 5.2 that dry olivine may sample an effective grain-boundary viscosity distribution spanning several orders of magnitude.

These results should be interpreted cautiously. The fit is not a unique inversion for grain-boundary properties, and it inherits the assumptions of the reduced-order model, including the adopted low-temperature elastic parameterization, the laboratory anchoring of the reference sliding time, the neglect of additional attenuation mechanisms, and the use of a single activation energy for all boundary types. Nevertheless, the exercise demonstrates that distributed-viscosity EAGBS is capable, in principle, of producing both a nonlinear high-temperature velocity reduction and attenuation levels of order $Q^{-1}\sim 10^{-2}$, making it a plausible contributor to upper-mantle seismic anelasticity.

Other candidate mechanisms have also been proposed to explain the missing dispersion and attenuation in the upper asthenosphere. One interpretation invokes poroelastic effects [mavkoVelocityAttenuationPartially1980, takeiEffectPoreGeometry2002] and enhanced anelasticity [faulShearWaveAttenuation2004] associated with partial melting. However, geochemical constraints suggest that the melt fraction in the upper mantle is very small [mckenzieConstraintsMeltGeneration2000], making it difficult to attribute the observed large dispersion to melt alone [takeiPhaseFieldModelingGrain2019]. Another possible explanation is grain-boundary premelting [yamauchiApplicationPremeltingModel2020a], motivated by the observation of a high-frequency attenuation peak in borneol rock-analogue experiments [takeiTemperatureGrainSize2014, yamauchiPolycrystalAnelasticityNearsolidus2016]. This peak exhibits scaling relationships different from those of the EAGBS peak considered in this study. However, the physical mechanism by which premelting enhances anelastic behaviour and produces such a peak remains unclear. Further theoretical and experimental work is therefore needed to clarify the role of premelt in the upper asthenosphere.