Hi, @lixinyv and @jbnaliboff,
Thank you very much for your insightful and helpful comments. Your suggestions guided my follow-up tests and clarified my confusion.
I have conducted further tests by varying the shear modulus (\mu). It indeed appears that as \mu increases, the numerical solution converges toward the viscous analytical solution, and the oscillations gradually subside (Fig. 1). This suggests that the oscillations are not accidental but are likely related to numerical instability in cases where \mu is relatively low.
Figure 1. Topographic evolution over time of different
\mu, based on the maximum topography value.
Figure 2. Variations of time, time step size (dt), and viscosity (viscoelastic viscosity) with simulation steps, related to Figure 1
Sorry for the confusion. This test case is actually adapted from an Underworld2 tutorial on viscoelastic half-space loading. I noticed a discrepancy: Underworld’s result appears smoother than ASPECT’s initial “drop” (Fig. 3, blue and green lines). Since I did not find any obvious initial stress supporting the topography in the Underworld example, now I think this difference might stem from how the two codes implement viscoelasticity, including the way they use the elastic time step size (dt_e, default 10).
When I reduced dt_e in Underworld to be closer to \Delta t_c (the actual time step, fixed at about 0.12), which enabled a more stringent capture of elastic effects, it also showed rapid subsidence and oscillations (Fig. 3, purple line). This behavior corresponds to the initial elastic effect observed in ASPECT. And this might indicate that ASPECT’s behavior reflects a more direct capture of the elastic response for the same dt_e condition.
Figure 3. Topographic evolution over time when
\mu=0.03 / 3 in ASPECT and Underworld, based on the maximum topography value.
A simplified viscoelastic analytical solution derived by AI (Gemini) is provided below (for reference only).
w(t) = w_m \left( \frac{1}{1 + \frac{\rho g \lambda}{4\pi \mu}} \right) \exp \left[ -t \left( \frac{\eta}{\mu} + \frac{4\pi \eta}{\rho g \lambda} \right)^{-1} \right]
[Correction]: The term \frac{\rho g \lambda}{4\pi \mu}, which represents the stiffness ratio, was inverted previously and has now been corrected.
As shown in Fig. 4, this solution matches ASPECT’s results with high precision across various \mu values, with minor deviations in the case of \mu=0.03.
Figure 4. Topographic evolution over time for different
\mu and related analytical solutions.
In summary, I agree with the view that these oscillations are part of the adjustment process and are to be expected. I would like to thank you again for your attention and guidance regarding this issue — your help has deepened my understanding of the viscoelastic property and its implementation in ASPECT. It is worth noting that, due to certain reasons, I have not yet calculated the magnitude of the dimensional shear modulus; however, the observed trend in \mu might point to this being a relatively extreme scenario. I’d be glad to participate if there are any further tests.
Best regards,
Junru
Analytical Derivation of Viscoelastic Half-Space Loading
Disclaimer: This derivation was generated with the assistance of an AI (Gemini) to serve as a qualitative reference for the discussion.
1. Problem Definition
We consider a semi-infinite half-space (y > 0) composed of a Maxwell viscoelastic material with density \rho, viscosity \eta, and shear modulus \mu. The surface is subjected to an initial periodic displacement w(x, 0) with no initial stress support (instantaneous loading):
w(x, 0) = w_m \cos(kx)
where k = 2\pi/\lambda is the wavenumber. We seek the time evolution of the surface topography w(x, t).
2. Constitutive Law (Maxwell Model)
For a Maxwell body, the total strain rate tensor \dot{\epsilon}_{ij} is the sum of the elastic and viscous strain rates:
\dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^e + \dot{\epsilon}_{ij}^v = \frac{\dot{\sigma}_{ij}}{2\mu} +\frac{\sigma_{ij}}{2\eta}
Applying the Laplace transform (f(t) \to \tilde{f}(s)) to the constitutive equation, where \partial/\partial t \to s:
s \tilde{\epsilon}_{ij} = \frac{s \tilde{\sigma}_{ij}}{2\mu} + \frac{\tilde{\sigma}_{ij}}{2\eta} = \tilde{\sigma}_{ij} \left( \frac{s}{2\mu} + \frac{1}{2\eta} \right)
This allows us to define an effective viscosity \eta(s) in the Laplace domain that relates stress and strain rate formally as \tilde{\sigma}_{ij} = 2\eta(s) \tilde{\dot{\epsilon}}_{ij}:
\frac{1}{\eta(s)} = \frac{s}{\mu} + \frac{1}{\eta} \implies \eta(s) = \frac{\mu}{s + \frac{\mu}{\eta}}
3. Application of the Correspondence Principle
According to the Correspondence Principle, the solution to a linear viscoelastic problem can be obtained from the corresponding elastic or viscous solution by replacing the constant material parameters with their s-dependent counterparts in the Laplace domain.
The characteristic equation for the relaxation rate s of a purely viscous half-space under gravity is :
\rho g + 2k \eta s = 0
(Note: Here 2k\eta s represents the viscous stress term balancing the gravitational restoring force \rho g.)
Substituting the viscous viscosity \eta with the Maxwell effective viscosity \eta(s):
\rho g + 2k s \cdot \eta(s) = 0
\rho g + 2k s \left( \frac{\mu}{s + \frac{\mu}{\eta}} \right) = 0
4. Derivation of the Relaxation Rate
Rearranging the characteristic equation to solve for s:
\rho g \left( s + \frac{\mu}{\eta} \right) + 2k \mu s = 0
s (\rho g + 2k\mu) + \frac{\rho g \mu}{\eta} = 0
The root s (which corresponds to the inverse relaxation time -1/\tau_{ve}) is:
s = - \frac{\frac{\rho g \mu}{\eta}}{\rho g + 2k\mu}
Dividing numerator and denominator by \rho g \mu:
s = - \frac{\frac{1}{\eta}}{\frac{1}{\mu} + \frac{2k}{\rho g}} = - \frac{1}{\frac{\eta}{\mu} + \frac{2k\eta}{\rho g}}
Thus, the effective viscoelastic relaxation time \tau_{ve} is:
\tau_{ve} = \frac{\eta}{\mu} + \frac{2k\eta}{\rho g} = \tau_{Maxwell} + \tau_{viscous}
Substituting k = 2\pi/\lambda:
\tau_{ve} = \frac{\eta}{\mu} + \frac{4\pi \eta}{\rho g \lambda}
5. Initial Elastic Response (t=0^+)
In the limit of s \to \infty (equivalent to t \to 0), the material behaves purely elastically. The initial surface displacement w(0^+) is less than the applied load shape w_m because the elastic stiffness supports part of the load instantaneously.
Using the balance of forces at t=0: F_{gravity} + F_{elastic} = 0.
The effective stiffness ratio implies:
w(0^+) = w_m \frac{2k\mu}{\rho g + 2k\mu}= w_m \left( \frac{1}{1 + \frac{\rho g \lambda}{4\pi \mu}} \right)
6. Final Solution
Transforming back to the time domain, the topography evolution is:
w(t) = w(0^+) \exp\left( -\frac{t}{\tau_{ve}} \right)
w(t) = w_m \left( \frac{1}{1 + \frac{\rho g \lambda}{4\pi \mu}} \right) \exp \left[ -t \left( \frac{\eta}{\mu} + \frac{4\pi \eta}{\rho g \lambda} \right)^{-1} \right]
