Viscoelastic half-space loading#
Benchmark of infinite viscoelastic half-space loaded/unloaded by an axisymmetric cylinder, using the solution of: Nakiboglu, S.M. and Lambeck, K., A study of the Earth’s response to surface loading with application to Lake Bonneville. Geophys. J. R. astr. Soc. (1982) 70, 577-620.
The benchmark compares their solution for the surface displacement following the application of a prescribed loading function (Eq. 27) with the deformation of the free surface in ASPECT after imposing surface boundary tractions.
A surface pressure of rho_lgH0 (where rho_l is the load density, H0 is the load height) is applied on the surface for r<r0 (where r0 is the load radius), for t>0. The load is either instantaneously loaded/emplaced at times t0<t<t1, or thinned linearly over this time interval (“linear unloading”). The load is fully removed by t=t1, and the surface relaxes. This is done both in a 2-D and 3-D geometry (by symmetry the load is centered on the left boundary or left/front corner). The input files are: ‘free_surface_VE_cylinder_2D_loading.prm’ (instantaneous loading) ‘free_surface_VE_cylinder_3D_loading.prm’ ( “ “ “ ) ‘free_surface_VE_cylinder_2D_loading_unloading.prm’ (linear unloading) ‘free_surface_VE_cylinder_3D_loading_unloading.prm’ ( “ “ )
Changing the stress instantaneously as the load is applied/removed leads to oscillations in the elastic response. This can be addressed by setting a fixed elastic timestep which is several times longer than the numerical (“Maximum”) timestep and turning on stress averaging: ‘free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm’ ‘free_surface_VE_cylinder_3D_loading_fixed_elastic_dt.prm’
The ‘topography’ output files may be compared against the analytical solution by running the script ‘compare_VE_soln.py’, which outputs the Nakiboglu and Lambeck (1982) solution for a cylindrical load on an infinite half-space. These outputs (‘soln…txt’) may be compared against the ASPECT ‘topography’ output files, using the provided gnuplot script (‘compare_viscous_def.gnuplot’ for maximum surface deflection through time, ‘compare_viscous_def_profile.gnuplot’ for deflection of profile through time). Equivalent scripts for the unloading case are also provided.
Note that while the analytical and numerical results for the deflection of the surface agree well near the center of the load (left boundary), the solutions do not match as well on the right (free-slip) boundary. The numerical free surface rides up vertically against the right boundary to conserve mass, while the analytical solution assumes an infinite half-space, predicting near-zero displacement far from the load center. This is also true for the viscous benchmark. This problem is less pronounced in 3-D as the extra mass may be distributed over a larger area. An open far (right/back) boundary resolves this problem in 3-D.
The solutions match well in 3-D. In 2-D, the geometries of the loading function are different (Cartesian in ASPECT vs cylindrical analytically). As such, the agreement in 2-D breaks down for small r0 (load width) or if the right boundary is placed further away.
The fit between the numerical and the analytical solution can be improved by choosing short elastic (and therewith even shorter numerical) time steps. The shorter the elastic time step is, the more instantaneous is the surface response to loading and unloading. When using a very small time step in the analytical solution it can be seen that the response is indeed instantaneous, favoring small elastic and computational time steps for an accurate numerical solution as well.