DOI: 10.31038/GEMS.2025773

Abstract

Laboratory experiments with rock samples show that creep at small strains is transient and is described by the linear hereditary rheological model of Andrade. Flows that restore isostasy (in particular, postglacial ones) cause deformations in the lithosphere that do not exceed 10^-3 and, therefore, demonstrate transient creep. The effective viscosity characterizing transient creep is lower than the effective viscosity at steady-state creep and depends on the characteristic time of the process under consideration. The characteristic duration of isostatic equilibrium recovery after the initial disturbance of the Earth’s surface relief does not exceed ten thousand years, and therefore the depth distribution of rheological properties differs from the distribution that corresponds to slow geological processes. The perturbations of the Earth’s surface relief caused by an initial small-scale perturbation that disrupts isostasy are considered. The propagation of vertical displacements along the Earth’s surface from the area of the initial perturbation is carried out by diffusion-type waves that arise during the process of isostasy recovery, and by convective waves that arise due to the vertical temperature gradient in the lithosphere. The solutions of the equations of continuum mechanics are obtained using the Fourier transform in the spatial horizontal coordinate and the Laplace transform in time.

Keywords

Transient creep, Isostatic recovery, Vertical movements of the Earth’s surface, Thermoconvective waves, Diffusion-type waves

Introduction

Andrade proposed a law to describe the transient creep he observed experimentally in 1910. This law was later confirmed in numerous laboratory studies of rock creep conducted at high temperatures and pressures typical of the Earth’s interior [1-5]. The well-known concept in mechanics that creep is transient at small deformations was first introduced into geophysics in the work of [6], where the idea was put forward that flows in the mantle associated with small deformations, and, in particular, postglacial flows, occur in the transient creep regime. The concept of transient creep at small deformations of geomaterial was further developed in the monograph [7]. Experimental and theoretical justifications for the applicability of the Andrade rheological model in studying dynamic processes in the lithosphere and mantle are presented in [8]. The lithospheric plate is a cold boundary layer formed by mantle convection, and the thickness of continental plates beneath cratons can exceed 200 km. The rheological model of a power-law non-Newtonian fluid, which describes steady-state creep and is commonly used in modern geophysical studies, leads to a very high effective viscosity characterizing creep at small deformations. With such a high effective viscosity, the lithosphere would be purely elastic even over geological times. Transient creep corresponds to a much lower effective viscosity than steady-state creep. Therefore, transient creep must be taken into account when examining geophysical processes in the lithosphere. The effective viscosity corresponding to transient creep depends on the characteristic duration of the geophysical process under consideration. In [8-11], low-amplitude thermoconvective waves in the lithosphere, the existence of which is due to transient creep, were considered, and it was shown that thermoconvective oscillations (standing waves) lead to the formation of sedimentary basins on continental cratons. The characteristic time of this process is about 108 years. The work [12] considers the recovery of isostatic equilibrium after an initial small – scale disturbance of the earth’s surface relief. As a result of the recovery process, the Earth’s surface returns to a flat position, which corresponds to the equilibrium state with a uniform horizontal density distribution. The characteristic duration of this process does not exceed 1000 years, and therefore the distribution of rheological properties by the depth of the lithosphere and crust differs from the distribution that corresponds to slower processes associated with convective motion. A process with a characteristic time of about 1000 years is a fairly fast process, the study of which does not require taking into account the influence of the vertical temperature gradient present in the lithosphere, since thermal effects are associated with very slowly occurring thermal conductivity. However, this process is slow enough to neglect the elasticity, compressibility, and inertia of the lithosphere.

In this paper, we will consider the movements of the Earth’s surface caused by initial vertical displacements. If a vertical displacement disrupting the isostatic equilibrium of the crust has occurred in some area of the Earth’s surface, then the process of isostasy recovery is reduced not only to a gradual decrease in the initial vertical displacements in this area, but also to the propagation of displacements beyond its boundaries. The processes of displacement propagation in the crust can be called diffusion-type waves. Waves of this type appear when solving parabolic equations of the theory of diffusion or heat conduction. These waves describe the process of propagation of temperature disturbance from the area where the disturbance occurred to the environment. The temperature gradually increases in remote areas, and the temperature in the initially disturbed area decreases. The propagation of the disturbance is accompanied by its strong attenuation, and over time, the temperature disturbance disappears everywhere. This is exactly how the process of displacement propagation occurs in the elastic crust, where the displacement acts as an analogue of temperature disturbance. It will be shown what movements of the Earth’s surface created by diffusion-type and thermoconvective waves caused by the initial vertical displacement of this surface.

Rheological Model

Laboratory studies show that at small deformations transient creep occurs, in which creep deformations depend linearly on the applied constant stresses

where  f(t) is the creep function, which provides an analytical description of transient creep, ε_ij and  is the strain tensor measured from the state at the moment of stress application. For rocks, the creep function at high temperatures is well described by Andrade’s law

where A is the Andrade rheological parameter. At short times, transient creep obeys Lomnitz’s law, but already at times of the order of a day, Andrade’s law becomes valid [13]. Therefore, Andrade’s law will be used in the study of low-frequency waves. The attenuation of high-frequency seismic waves is described by Lomnitz’s rheology [14], but in this paper, although we will talk about Rayleigh surface elastic waves, their attenuation will not be considered. To generalize the results of experiments carried out at constant stresses to the case of variable stresses, we can use the linear Boltzmann theory, valid for sufficiently small deformations. This theory leads to an integral relationship between strains and stresses

where t is the observation time, and K(t) is the integral creep kernel determined by the creep function

As follows from (2) and (4), the creep kernel corresponding to Andrade’s law has the form

The rheological model described by equations (3) and (5) will be called the Andrade model. This model generalizes Andrade’s law to the case of variable stresses. At sufficiently large deformations, transient creep is replaced by steady-state creep, which is described by the rheological model of a power-law non-Newtonian fluid. The value of the Andrade rheological parameter depends on temperature, pressure, and mineralogical composition. The depth distribution of the Andrade parameter was obtained in [11]. Since the Andrade parameter decreases with increasing temperature, and the temperature increases with depth in the lithosphere, this parameter decreases with depth. In the upper crust, the thickness of which is about 20 km, the value of this parameter is estimated as A ≈ 10¹⁶ Pa∗ s^1/3 . With such a high value of  A, the upper crust behaves as an elastic layer even on times of the order of the age of the Earth. More precisely, the upper crust, the thickness of which is about 20 km, is brittle-elastic, and elasticity dominates only in the lower layer of the upper crust, the thickness of which is about 10 km. In the uppermost layer of the crust, brittleness dominates, and its strength is very low [11]. We will assume that the upper crust with a thickness of 20 km is an elastic layer with a reduced shear modulus due to brittleness. Within the framework of the simplified model of the lithosphere beneath the upper crust used in the paper,we use the depth-averaged estimate A ≈ 10¹² Pa∗ s^1/3.

Since the total deviatoric deformation of the medium can be represented as the sum of the deviatoric creep deformation (1) and the deviatoric elastic deformation

where μ is the elastic shear modulus, an elastic-creeping medium whose creep is transient is described by the equation

In linear stability theory, it is assumed that the strains and stresses depend on time as exp(λt), where λ is the complex decrement, and the right-hand side of equation (3) takes the form

where the asterisk denotes the Laplace transform, which is used here only to evaluate the integral in equation (8). The Laplace transform of the creep kernel (5) yields

where the gamma function is Γ (1/3)≈3. Linear stability theory considers the behavior of a mechanical system at large times elapsed since the occurrence of a small disturbance. Therefore, in equation (8), the upper limit of integration is t=∞.

Thus, when the time dependence has the form exp(λt), at large times the effective shear modulus of the Andrade medium has the form

and the effective Newtonian viscosity is written as

where  is the characteristic time of the process. As follows from (11), on time scales of the order of 1000 years, small-scale postglacial flows, which are characterized by an average value of the rheological parameter A ≈ 10¹² Pa∗ s^1/3, correspond to an effective viscosity ℑ_eff ≈ 10¹⁹ Pa∗ s. Such estimate agrees in order of magnitude with the viscosity estimates found when considering small – scale postglacial flows within the framework of the rheological model of a Newtonian fluid [15]. Consequently, the estimates of Andrade parameter obtained on the basis of data from laboratory experiments with rock samples at temperatures and pressures characteristic of the Earth’s interior [11], correspond to the observational data used to estimate the viscosity of the upper mantle [15]. The effective viscosity for the Andrade medium depends on the characteristic time of the process under consideration and, therefore, the effective viscosity found for postglacial flows with a characteristic time of 1000 years cannot be used in the study of slower geological processes [8,13].

Thermoconvective and Elastic Surface Waves

We consider an elastic thin plate lying on a layer with Andrade rheology. The origin of coordinate is placed on the upper surface, and the z-axis is directed vertically upward. The thin plate (z = 0) models the upper elastic crust, and the layer (-1 < z < 0) – the underlying lithosphere. The equations describing the disturbances of the lithostatic equilibrium of an incompressible medium are written as

where p  is the pressure perturbation, σ_xx, σ_xz and σ_zz are the components of the deviatoric stress tensor,  v_x and v_z are the velocities, θ is the temperature perturbation,  x is the thermal diffusivity,  is the vertical gradient of the unperturbed temperature, which is assumed to be uniform throughout the depth of the lithosphere, ρ is the density, α is the coefficient of thermal expansion, g is the gravitational acceleration. The velocities, stresses, pressure and temperature perturbations are functions of the vertical spatial coordinate, the horizontal coordinate x, and time t. Equations (12) and (13) describe the two-dimensional motion of the medium (the motion occurs in the xz plane) taking into account the Archimedes force. Equation (14) is the condition of incompressibility of the medium, and equation (15) is the heat balance equation. Equations (12) – (15) use the Boussinesq approximation, within which the mechanical compressibility of the medium can be neglected, and thermal compressibility can be taken into account only in the equations of motion. Rheological equations are added to equations (12) – (15)

where F(λ) is the complex viscosity of the Andrade viscoelastic medium (7)

Equations (16) relate deviatoric stresses to strain rates, which are defined as

On the upper surface of the lithosphere (z = 0) the boundary conditions, determined by the force action of the elastic plate, are imposed

where ν is Poisson’s ratio, N is the flexural rigidity of the elastic plate with thickness h. The displacements  u_x and u_z in the plate are equal to the displacements in the underlying layer at z = 0. The temperature boundary condition is added to the boundary conditions (19) and (20)

The incompressibility condition for a viscous medium, under which equation (14) is valid, is written as

where τ is the characteristic time of the flow under consideration, and K is the bulk modulus. It follows from (23) and (11)

In equations (12) and (13), the inertial terms can be neglected under the condition

from which it follows

Slow flows with a large characteristic time τ are called creeping. Convective flows in the Earth are a typical example of creeping flows. The lithosphere, which exhibits both elastic and viscous properties, is described by the Maxwell viscoelastic rheological model. Elasticity can be neglected for a sufficiently slow flow

where µ is the elastic shear modulus, and η_eff /μ is called the Maxwell time. It follows from condition (27)

To move to dimensionless variables, we introduce the following scales: the length scale is the layer thickness d, the time scale is d²/x , where κ is the thermal diffusivity, the velocity scale is κ/d, the pressure (and stress) scale is kη_A /d², where η_A is the viscosity scale for the Andrade medium. The temperature difference between the hot lower and cold upper surfaces of the layer is taken as the temperature scale. For the Andrade medium η_A = A(d² /x)^2/3, and then with an exponential dependence on time, the dimensionless effective viscosity takes the form

and the Rayleigh number for the Andrade medium is defined as

The lithosphere is characterized by the following values of physical parameters [16]:

The thickness of the continental lithosphere, which is taken as the length scale, is estimated as d = 2∗ 10⁵ m. Since the Andrade parameter for the lithosphere is estimated as 10¹² Pa∗.s^1/3 m, the Rayleigh number, according to (31), is estimated as Ra ≈ 80.

We will represent the vertical velocity as

where λ is the complex increment, k is the real wave number. We will represent all the physical variables in a similar form. This representation allows us to reduce the system of partial differential equations (12) – (18) to a system of ordinary differential equations, in which all variables characterizing the strain rates, stresses, and pressure depend only on the vertical coordinate z. The characteristic time τ for a flow with a decrement λ can be represented as .

Passing to dimensionless variables and substituting relations (32) into equations (12) – (15), (16) and (18), we obtain relations connecting the amplitudes of pressure, temperature, and components of the deviatoric stress tensor with the amplitude of the vertical velocity

Excluding from the equations the amplitudes of all physical variables, except for the amplitude of the vertical velocity V_z , we arrive at an ordinary differential equation, valid for small values of λ, for which we can neglect elasticity and inertia of the layer (-1 < z <0) simulating the lithosphere,

From (34) it follows that the influence of the temperature gradient can be neglected when

The boundary conditions on the upper deformable surface of layer z = 1 are

where U_z is the vertical displacement of the upper boundary of the layer, i.e. the deviation of the boundary surface from the plane z = 1. Equations (36) and (37) follow from the condition of equilibrium of the forces acting on a unit area of the disturbed surface of the layer, and equation (38) follows from the condition of vanishing of the temperature disturbance on this surface. Equations (36) and (38) include the displacement of the boundary U_z due to the fact that in the state of equilibrium there is vertical gradients of pressure and temperature in the layer. Condition (38) assumes that the boundary motion is determined by the motion of material points located on this boundary. In equations (36) and (37), dimensionless parameters are introduced:

which characterize the deformable surface. For sufficiently small k, the influence of a thin elastic plate modeling the upper crust is negligible.

Using the parameter φ, the Rayleigh number can be written as  The parameter φ describes the mobility of the boundary and is the ratio of the additional hydrostatic pressure caused by the boundary displacement to the characteristic viscous stress in the layer. The smaller the parameter φ, the more mobile the boundary. For a very large value of φ, the upper surface of the layer behaves like a fixed boundary.

Under each lithospheric plate, there is an isothermal core of large-scale mantle convection, and there is movement with a constant horizontal velocity at the lithosphere-mantle boundary [16]. Therefore, in the coordinate system that moves together with the lithospheric plate, the following conditions must be imposed at the lower boundary of the plate

Thus, the lower boundary of the lithosphere, considered as a boundary layer of large-scale convection, is isothermal (the temperature perturbation is zero at this boundary) and “solid” (zero velocity components at this boundary).

The general solution of the ordinary differential equation (34) is written as

where   q₁ , q₂ , q₃ are the roots of the equation

Investigating surface waves, for which disturbances do not penetrate into the deepest layers of the lithosphere, we consider only sufficiently large values of the wave number (k > 3), for which  e^{q1 z} , e^{q2 z} and e^{q3 z} are close to zero at the lower boundary (z = -1). In this case, to satisfy the boundary conditions at the lower boundary of the lithosphere, it is sufficient to set C₄ = C₅ = C₆ = 0 Substituting (40) into the boundary conditions (36) – (38) on the upper surface, we arrive at a system of three homogeneous equations for three unknowns C_{1 ,} C₂ and C₃. In order for this system to have a solution, it is necessary to equate its determinant to zero

As a result, we obtain the characteristic equation relating the decrement λ to the wave number k. For large wave numbers (k > 5), the imaginary parts of the complex decrements are very small. This means that the temperature gradient present in the lithosphere has a very weak effect and the Rayleigh number in equations (34) and (41) can be neglected. For k = 3, the decrement takes the value . Such a convective wave propagates with the dimensionless velocity . The dimensional velocity of this wave is extremely small: m/year. Values of k < 3 cannot be considered since equation (42) was obtained for surface waves. However, the study of thermoconvective waves [8,9,10,11] allows us to assume that for k < 3 the imaginary part of the decrement and the propagation velocity increase, while the real part of the decrement (attenuation) tends to zero.

The characteristic equation (42) is obtained for sufficiently small values of k and λ, for which the elasticity and inertia of the medium can be neglected. For large values of the wave number k and the decrement λ, the viscosity of the medium and temperature effects can be neglected but the elasticity and inertia of the medium must be taken into account. In this case, it is more convenient, while maintaining the length scale d=2∗ 10⁵ m, to move to the velocity scale  and the time scale Then λ=iw (the decrement becomes purely imaginary), and the characteristic equation takes the form

This equation has solutions w(k) = Vk, where V ≈ 0.95529 is the dimensionless velocity of the surface wave. Its dimensional velocity is slightly lower than the velocity of the bulk transverse wave. Taking into account the effect of gravity (non-zero parameter φ) slightly increases the frequencies ω(k) ≈ 0.95531k. Equation (43) describes the Rayleigh wave in an incompressible medium. Taking into account compressibility, equation (43) becomes

where  This parameter can be considered small, since it has little effect on the result. The frequencies ω found from equation (43) and from equation (44) differ by approximately 3%.

When k and λ are small enough to neglect elasticity and inertia, but not small enough to take into account thermal effects, the complex decrement λ turns out to be a real number (lmλ = 0). In this case, which will be discussed in the next section, there is no wave motion at a fixed k but we can speak of a diffusion – type wave in the presence of a whole spectrum of values of k.

Laplace Transform. Diffusion-type Wave

Solving the problem of excitation of surface waves by the initial disturbance of the Earth’s surface, we will use the Laplace transform with respect to time t and the Fourier transform with respect to the spatial horizontal coordinate x. Representation (32) introduces the Fourier transform but not the Laplace transform, which allows us to find not only solutions of the form exp (λt).

The Laplace transform is a convenient mathematical apparatus for solving differential equations with initial conditions. The Laplace transform f∗(s), denoted by an asterisk, is related to the original f(t) by the equation

The Laplace variable s is a complex number, unlike the Fourier variable k, which is a real wave number.

Using the Laplace transform, we can write the rheological equation (3) as

where G∗_A  is the Laplace analogue of the shear modulus for the Andrade medium, Γ(1/3) ≈ 3 is the gamma function. Equation (7) corresponds to the Laplace image

where G∗ is the Laplace analogue of the shear modulus for an elastic-creeping medium.

As follows from (46), the elasticity of the medium can be neglected if the condition

is satisfied, under which the Laplace analogue of the shear modulus for the Andrade medium G∗_A is significantly less than the elastic shear modulus μ. According to the property of the Laplace transform, it follows from (47)

The compressibility of the medium can be neglected if the Laplace analogue of the shear modulus for the Andrade medium G*_A is significantly less than the elastic bulk modulus K. Since the elastic shear modulus μ is less than the bulk modulus (K ≈ 3μ), condition (48) allows us to neglect not only the elasticity but also the compressibility of the medium. Since μ ≈ 6. 10¹⁰ Pa and A ≈ 10¹² Pa ∗ s^1/3 , it follows from (48) that the lithosphere beneath the elastic upper crust behaves as a creeping Andrade medium, not exhibiting elasticity or compressibility at times exceeding 10⁴ s.

The effect of inertia is negligibly small (inertial forces are small compared to the forces arising during deformations of the Andrade medium), if

The right-hand side of (49) is estimated as 20 s for the lithosphere. Consequently, neglecting elasticity, inertia can be neglected even more so. Isostatic recovery processes with characteristic times not exceeding 10,000 years are slow enough to neglect elasticity, compressibility, and inertia, but are not slow enough to take into account the buoyant Archimedean force, the presence of which is due to the vertical temperature gradient in the lithosphere. The physical variables in equations (12) – (14) depend on the horizontal coordinate x, the vertical coordinate z, and time t. Applying the Fourier transform in the coordinate x and the Laplace transform in time t to these equations and excluding all physical variables except the vertical displacement, we arrive at the relations

where the differential operator D = d/dz is introduced, and U₀=U₀(x, z) is the initial (t = 0) distribution of vertical displacements. In the equations (50) – (53), the Fourier transforms of the physical variables are denoted by corresponding capital letters, the Laplace images are marked with the asterisk, k is the wave number (the Fourier variable), s is the Laplace variable. Then we obtain the ordinary differential equation for the vertical displacement

The solution of equation (54), which satisfies the boundary condition for z⟶ – ∞, is written as

where C₁ and C₂ are arbitrary integration constants that depend on the Laplace variable, and the wave number k can take negative values. Substitution of the solution (55) into the boundary conditions (19) – (20) allows us to eliminate arbitrary constants and represent the vertical displacement of the upper surface (z = 0) in the form

where G^*_A = As^1/3 is the analog of the shear modulus for the Andrade medium. It should be noted that the substitution of (50) and (53) into the boundary condition (19) leads to the relation C₂ = – C₁ /k, at which horizontal displacements and tangential stresses are absent on the upper surface. Thus, when searching for a solution to equation (54), one can replace the boundary condition (19) with the condition U^*_x= 0 at z = 0 or with the condition Σ^*_xz= 0 at z=0.

In order to find the asymptotic (small times) dependence of the Laplace origin on time, it is sufficient to expand the Laplace image in a series in powers of s in a neighborhood of s = ∞ and invert by Laplace each term of the series [17]. The right-hand side of (24) for s⟶ ∞, is representable in the form of a series

Inverting the terms of the series (57), we obtain the asymptotic dependence of the vertical displacement on time

where the gamma-function at the point 4/3 is estimated as

The asymptotic dependence (59) is valid when

As follows from (60), the asymptotic dependence (59) can be represented in the form

The function Φ(k) has a sharp minimum, which is found from the condition

If we switch to the length scale introduced above d=2∗10⁵ m, this wave number takes the value K_m= 3.6

At large times, the displacements of the surface are described by another asymptotic formula. In the neighborhood of the point s = 0, the right-hand side of equality (56) can be represented as a series  

According to the theorem on the asymptotic behavior of the original [17], the Laplace original at large times can be represented as a series whose terms are obtained as a result of the inverse Laplace transform of each term of the series (64). Keeping only the first term of the expansion, we find

Asymptotic dependence (65) is valid at large times, when

As follows from (65), at large times, harmonics with different wave numbers k decay according to the same law t^-1/3 and the effect of propagation of surface displacements, which occurs at small times, disappears. Since the minimum value of Φ(k) is achieved at k=k_m , inequality (66) is satisfied for any wave numbers if

If the medium underlying the elastic upper crust had the rheology of a Newtonian fluid with viscosity η, the analogue of the shear modulus G^*_A(s) should be replaced by ηs, and equation (56) would be written as

Inversion of the Laplace image (68) yields

where τ is the recovery time, depending on the wave number k. Relation (69) is valid for any times t, in contrast to relation (61), which characterizes the Andrade medium and is valid only for not too large times, limited by condition (60).

Let at the initial instant t = 0 the displacement of the upper surface (z = 0) be given as

The Fourier transform of the function (70) is

and the inverse Fourier transform gives

Inversing the Fourier image (61), we find

In the case when u₀=x is an even function, (73) can be written as

Let the initial displacement have the form of a “step”

If x<-1/2 or x > 1/2 

The Fourier transform for such a “step” has the form

For a sufficiently large width l of the initial perturbation, the values of the function (76) are very small when k > 2π /l, i.e., the wider the perturbation region, the narrower the range of wave numbers k, in which the Fourier image is different from zero. Thus, the integration on the right-hand side of (74) is carried out over the region k<2π/l in the neighborhood of the point k = 0. According to the solution (55), for a fixed wave number k, the isostatic flow causes displacements in the lithosphere, depending on the depth as exp (- kz). Since k < 2π /l, we can say that this flow penetrates into the lithosphere to a depth of the order of l/π. Therefore, the flows arising after the removal of small-scale glaciations or other surface loads (for example, drying salt lakes), for which l does not exceed 200 km, are concentrated in the lithosphere. Caused by large-scale glacial loads (l ≈ 1000 ÷ 3000 km) flows that penetrate into the low mantle and recover isostasy over a period of time of about 10,000 years, are not considered in this paper.

As follows from (76), when the width of the initial perturbation is small (lk << 1), the Fourier transform does not depend on k

The image (77) corresponds to a point initial perturbation

where δ (x) is the delta function

In the case of a point initial perturbation (perturbation of any initial width l can be regarded as a point perturbation when we consider displacements at a sufficient distance from the initial perturbation, that is, for x >> l) , it is possible to obtain an analytic solution of the problem of vertical surface motions. For the point initial perturbation, the dependence of the vertical displacements of the surface on the horizontal coordinate and time is determined by the integral

As follows from (62), the expansion of the function Φ(k) in a power series in the neighborhood of has the form k=k_m

where

The power series (80), which represents the function Φ(k) given by (58), converges when |k-k_m|≤k_m, i.e., the radius of convergence of this series is R=k_m

After changing the variable

the integral (79) takes the form

Equation (83) can be rewritten as

where

The integral on the right-hand side of equality (84) is calculated using the saddle point method used in the theory of functions of a complex variable [18]). The stationary point v_o is found from condition

The function of the complex variable f(v), which is given by (85), has a stationary point

As follows from (87), in a neighborhood of the stationary point v_o, the function f (v) can be represented in the form

Substituting (88) into (84), we obtain

By choosing such a path of integration in the complex plane, which is determined by the saddle-point method, we find the relation

where erf(x) is the error function, and R=k_m is the radius of convergence of the power series (80). At times significantly exceeding 10¹⁰s ≈ 300 years, the condition

It is known that erf(x)=1 for x>>1, and the value erf(x) is close to 1 for x > 1.

Thus, as follows from (79), (89) and (90), the required distribution of the vertical displacements of the surface takes the form

The found solution (92) is valid for sufficiently long times (from several hundreds to several thousand years). The upper bound on time is imposed by condition (60), in which k=k_m.

The graphs in Figure 1 are constructed by the relation (92) and show the dependence of the vertical displacements on the horizontal coordinate at different instants of time.

Figure 1: The dependence of the vertical displacements of the Earth’s surface on the horizontal coordinate at different times for the case when the width of the region of the initial displacement l is 10 km. Curve 1 corresponds to 30 years, curve 2 to 300 years, curve 3 to 1000 years.

By differentiating the right-hand side of (92) with respect to t for a fixed value of x, it is not difficult to find the velocity of the vertical motion of the Earth’s surface at points sufficiently far from the region of the initial disturbance of the relief. For example, when u_o= 100 m, l = 10 km, x = 100 km, the velocity  reaches its maximum value (about 1 mm/year) during the time t ≈ 600 years. The generalization of the considered case of the initial point vertical displacement to the 3D problem formulation implies not too much new to understanding of the process under consideration: it is sufficient to replace the coordinate x by the polar coordinate r in Figure 1 (the solution does not depend on the polar angle for the point perturbation of vertical displacement). The initial disturbance of the Earth’s surface leads not only to diffusion-type waves, but also to seismic waves and convective waves. Diffusion-type waves cannot be considered without the initial disturbance, which produces a spectrum of harmonics. Each of these harmonics is not a traveling wave. The wave effect appears only when the entire spectrum is considered.

In the previous section, traveling harmonic waves, both convective and elastic, were considered. The solutions obtained can be represented as

where u_o is the initial vertical displacement of the surface, ω=Imλ, Λ=Reλ, and λ is the found value of the complex decrement (Reλ<0). The use of the Laplace transform allows us to find an arbitrary constant C. To do this, we transform the original equations using the Laplace transform, and substitute the resulting general solution, containing arbitrary constants, into the boundary conditions Laplace transformed. As a result, we find the function F(s, k) and obtain a solution to the problem in the form of a Laplace image

To invert the Laplace image (94), we can use the well-known theorem on the asymptotic behavior of the original [17]. According to this theorem, in order to find an asymptotic solution at large times (t → ∞), it is sufficient to know the Laplace image in the neighborhood of the singular point . In our case, this singular point is a first-order pole, in the neighborhood of which

and the Laplace original has the form

The pole s_o is equal to the complex decrement λ found in the previous section.

Expression (96) is the Fourier image. To pass to the Fourier original, it is necessary to specify the function U_o(k) determined by the initial disturbance. As a result, we can obtain a solution in the form of a running wave packet.

Conclusion

The paper considers disturbances of the Earth’s surface relief caused by an initial small-scale disturbance that disrupts isostasy. The solutions of the equations of continuum mechanics are obtained using the Fourier transform with respect to the spatial horizontal coordinate and the Laplace transform with respect to time [19]. The propagation of vertical displacements along the Earth’s surface from the region of the initial perturbation is carried out by surface waves. At short times, the presence of vertical temperature gradient in the lithosphere can be neglected and a solution in the form of a decaying seismic wave can be obtained. In this case, the creep of the medium, leading to attenuation, is described by the rheological Lomnitz’s law. At very short times, creep can also be neglected, obtaining a solution in the form of an elastic Rayleigh wave. At very long times, the elasticity and inertia of the lithosphere can be neglected but, taking into account the vertical temperature gradient, a solution in the form of a thermoconvective wave can be obtained. At not too large and not too small times, one can neglect not only elasticity and inertia, as in the description of a thermoconvective wave, but also thermal effects, as in the description of Rayleigh waves, obtaining a wave of the diffusion type that occurs in the process of isostasy recovery.

The initial disturbance determines the spectrum of wave numbers k. Elastic Rayleigh waves are characterized by large values of k and propagate very quickly, determining the relief of the Earth’s surface only in the first seconds after the disturbance occurs. Thermoconvective waves are characterized by small values of k (large wavelengths) and propagate very Diffusion-type waves correspond to values of k lying between the wave numbers characteristic of thermoconvective waves and the wave numbers characteristic of Rayleigh waves. The wave character of the solution at short times is due to inertia, and at long times – to the vertical temperature gradient. In the case of Rayleigh and thermoconvective waves, each harmonic moves according to the law , where ω is the frequency and Λ(k) is the decrement. In a diffusion-type wave, each harmonic is motionless, and the harmonic with the wave number (wavelength decays most slowly. Wave motion appears after the inverse Fourier transform due to the dependence of attenuation on k.

Conflict of Interest

The author declares no conflicts of interest in this paper.

References

1. Goetze C (1971) High temperature rheology of Westerly J Geophys Res 76: 1223-1230.
2. Goetze C, Brace WF (1972) Laboratory observations of high-temperature rheology of rocks. Tectonophysics 13: 583-600.
3. Murrell SAF, Chakravarty S (1973) Some new rheological experiments on igneous rock at temperatures up to 1120°C. Geophys J Roy Astron Soc 34: 211-250.
4. Murrell SAF (1976) Rheology of the lithosphere – experimental Tectonophysics 36: 5-24.
5. Berckhemer H, Auer F, Drisler J (1979) High-temperature anelasticity and elasticity of mantle peridotite. Phys Earth planet Inter 20: 48-59.
6. Weertman J (1978) Creep laws for the mantle of the Phil Trans Roy Soc London A 288: 9-26.
7. Karato S (2008) Deformation of Earth An Introduction to the Rheology of Solid Earth. Cambridge university press.
8. Birger BI (1998) Rheology of the Earth and thermoconvective mechanism for sedimentary basins formation. Geophys J Inter 134: 1-12.
9. Birger BI (2000) Excitation of thermoconvective waves in the continental Geophys J Inter 140: 24-36.
10. Birger BI (2012) Transient creep and convective instability of the Geophys J Inter 191: 909 – 922.
11. Birger BI (2013) Temperature-dependent transient creep and dynamics of cratonic Geophys J Inter 195: 695-705.
12. Birger BI (2018) Andrade Creep at the Isostasy Recovering Flows in the Izv Phys Solid Earth 54: 849-858.
13. Birger BI (2007) Attenuation of seismic waves and universal rheological model of the Earth’s Izv Phys Solid Earth 43: 635-641.
14. Birger BI Internal friction in the Earth’ crust and transverse seismic waves. AIMS Geosciences 8: 84 – 97.
15. Cathles LM (1975) The viscosity of the Earth’s mantle. Princeton university press.
16. Turcotte and Schubert G (1982) Geodynamics: Application of Continuum Physics to Geological Problems. New York: John Wiley.
17. von Doetsch G (1974) Introduction to the Theory and Application of the Laplace Transformation.
18. Copson ET (2004) Asymptotic Cambridge University Press.
19. Birger I. (2024) Instability of the Earth’s lithosphere. Cambridge Scholars Publishing.

DOI: 10.31038/GEMS.2025772

Abstract

We show that the loss of land ice mass is determined not by global warming but by seismic activity, and thus neither supports nor rebuts global warming theories.

Keywords

Loss of land ice mass

There are lots of articles claiming that land ice mass loss is caused by global warming. Figures 1-3 show that the distribution of land ice loss in Greenland and Antarctica mimics the distribution of coastal seismic activity, strongly suggesting that the primary factor affecting land ice loss is not temperature but coastal seismic activity.

Figure 1: Map A shows GRACE and GRACE-FO observations of Greenland land ice mass change in 2002 – 2023, according to NASA. Map B shows magnitude ⩾ 4.0 earthquakes in 60.3°N-84°N, 70°W-10°W in 1997/1/1 – 2025/6/1. Map A is shown in the projection employed by NASA, while map B is shown in the projection employed by USGS. For the ease of comparison, each map shows the towns of Savissivik, Aasiatt, Nanortalik, Kulusuk, mount Gunnbjorn, points 77°N, 24°W, 82°N, 24°W, 82°N, 66°W, all marked with asterisks, as well as Scorseby Sound. Coastal earthquakes are marked with different colors, while quakes removed from the coast are shown in gray. The regions of large ice loss around Savissivik and Aasiatt and the region between the two towns with somewhat lesser ice loss correspond to earthquakes marked purple. The region of large ice loss around Nanotalik corresponds to earthquakes marked green. The region of large ice loss around Kulusuk, Gannbjorn corresponds to earthquakes marked blue. The region of medium ice loss between points 77°N, 24°W, 82°N, 64°W corresponds to earthquakes marked orange. The region between Scorse by Sound and 77°N, 24°W shows no ice loss, nor does it show any earthquakes. The region between points 82°N, 64°W and 79°N, 66°W shows only small ice loss, and only one earthquake. Map C is a copy of map be showing two largest landslides; the 2017/6/17 landslide was just next to two purple quakes, while the 2023/9/16 one was close to the two quakes marked brown.

Figure 2: MapA shows GRACE and GRACE-FO observations of land ice mass change in Antarctica, according to NASA. Map B shows magnitude ⩾ 4.0 earth-quakes south of 62°S in 1997/1/1-2025/6/1. Map A is presented in the projection selected by NASA, while map B is presented in the projection selected by USGS. Images C and D are parts of map A placed under the corresponding portions of map B. The largest rate of ice loss is on the Amundsen Sea sector. The second largest rate of ice loss is on the Antarctica Peninsula and Alexander Island 71°S, 70°W, with third largest rate of ice loss in Antarctica is in Queen Mary Land and Wilkes Land. The fourth largest rate of ice loss in Antarctica is Cape Andreyev. Map E shows magnitude ⩾4.0 earthquakes south of 62°S in 2015/1/1 – 2025/6/1. It shows no earthquakes in Queen Mary Land and Wilkes Land, the region has experienced most ice gain in 2021-2023. Map F shows all magnitude ⩾5.3 earthquakes in 1900 – 2024 south of 64°S; it reveals that the centers of coastal seismic activity practically coincide with the centers of ice loss.

Remarkable is the presence of antipodal symmetry in the phenomena discussed, As Figure 3 shows, the Arctic sinkholes of the past 30 years appeared almost antipodal to the centers of Antarctic ice loss. The largest loss of ice mass in Antarctica occurs in and close to Amundsen Sea sector, almost antipodal to Taymyr, Kara Sea, and Novaya Zemlya.

Figure 3: Map A shows the 2017/6/17 landslide in Nuugaatsiaq, 71.535°N 53.2125°W and the 2023/9/16 landslide in Dickson Fjord, 72.833°N, 26.95°W from Figure 1-C, marked by diamonds; recently-formed sinkholes in the Gyda, Yamal peninsulas, as well as one in Taymyr peninsula at ≈75.5°N, 108°E, marked by asterisks ⋆; an unusual melting in Auyuittuq National Park 67.883°N, 65.017°W in the summer of 2008, marked by a four-point star; the northernmost volcano Beerenberg; and the most powerful earthquake north of 64°N. Map B shows the southernmost volcano Erebus; three most powerful earthquakes south of 60°S; and the regions of ice loss from Figure 3. Map C shows the antipode of Antarctica contour superimposed on the contours of the Arctic. The two volcanoes are almost antipodal to each other, as are the 1933/11/20 and 1998/3/28 earthquakes. The Amundsen Sea sector of ice loss is almost antipodal to the recently-formed sinkholes Gyda and Yamal peninsulas; as well as the Novaya Zemlya and Severnaya Zemlya archipelagos, which, according to the University of Edinburgh, experienced the largest loss of ice in the Russian Arctic in 2010 – 2018. Antarctica Peninsula and Alexander Island 71°S, 70°W, which showed the second largest rate of ice loss, are almost antipodal to the recently-formed sinkhole in Taymyr. Queen Mary Land and Wilkes Land, which showed the third largest rate of ice loss, is almost antipodal to the 2017/6/17 landslide and Auyuittuq National Park. The fourth largest rate of ice loss in Antarctica is almost antipodal to the 2023/9/16 landslide and the giant oods around 71.08°N,26.83°W in the Scoresby Sound 70.5°N,25°W.

WL-QML, almost antipodal to Ban Island and Greenland. The Taymyr-Kara Sea-Novaya Zemlya region has been also marked by recently formed sinkholes. Figure 3 suggests that the much-larger- than-average loss of ice mass and the appearance of sinkholes are due to subglacial/subpermafrost thermal activity most likely caused by seismicity. Figures 1-3 confirm that the regions of ice loss mimic seismic activity in Antarctica as well as Greenland; the events in Greenland and Baffin Island marked in Figure 3A occurred near the regions of increased seismicity in Figure 1. NASA does not provide any information about ice mass loss along the arctic boundary of Russia, however, Figure 4 suggests that the recently-formed sinkholes are also related to seismic activity. That ice loss, in one form or another, may be caused by quakes, contemporaneous or precedent, is supported by Figure 5 and the 1958/7/10 UTC time (1958/7/9 local time) massive landslide caused by a magnitude 7.8 quake.

Figure 4: Map A shows magnitude ⩾ 3.9 earthquakes in 66°N−80°N, 34°E−180°Ein 2003/6/1 – 2025/6/1 along with sink-holes discovered after2015/5/1. Map B zooms in on the Taymyr peninsula. Map C shows all nuclear explosions in the region in 1973-1990, there have been no nuclear explosions in the region since. Map D is just a part of Figure 3C. The sink holes in Yamal are within the triangle formed by the 3 earthquakes around it, but the sinkhole in Gyda are not. However, the sinkholes in Gyda are just north of nuclear explosions shown. The Taymyr sinkhole, it seems to be a harbinger of the earthquakes to hit Taymyr shortly. Prior to 2003/6/1, only two quakes of that magnitude are known to have hit Taymyr, one on 1986/5/19, the other one on 1990/6/9.

Figure 5: On 2025/5/28, a huge portion of a glacier in the Swiss Alps had broken from the mountainside and crashed onto the village of Blatten at 46.417°N, 7.817°E, shown by an asterisk, approximately 16 km from the epicenter of the 1946/1/25 magnitude 6.2 quake at 46.499°N, 7.644°E, shown by a disk. Although the earthquake struck in 1946, smaller quakes in the region have not stopped until now; e.g. 2016/10/24 46.421°N 7.576°E magnitude 4.4, 2024/6/4 47.085°N 8.796°E magnitude 4.2, etc., not to mention numerous magnitude ⩽4.0 quakes.

DOI: 10.31038/GEMS.2025771

Abstract

During the Raman study of gray cast iron, we found graphite flakes decorated with spherical diamond crystals. The first-order Raman bands of these diamonds and the graphites that belong to them are given. The diamond pattern leads us to assume that the diamonds were primarily crystallographically arranged. The now arcuate arrangement lets us presume that stress deformed the graphite flake perpendicular to the c-axis. The often present Raman forbidden 867 cm^-1 line, which is usually IR-active, is found together with frequently found methane and benzene in graphite, maybe as intercalates. As a rule, the studied diamonds always show the graphite line besides the diamond band because the cross-section for sp₂ is much greater than for sp₃– bonded structures.

Keywords

Raman spectroscopy, Graphite flakes, Diamond in cast iron, Pressure questions, Forbidden Raman line, CH₄ intercalation, Deformation of graphite

Introduction

During the study of gray cast iron in “as-cast state”, we found diamond crystals in graphite for the first time using Raman spectroscopy. In the past, Sobolev et al. (1993) [1] synthesized diamond in cast iron by detonation. Forerunners were DeCarli PS, Jamieson JC (1961) [2] and Cowan et al. (1968) [3], who produced diamonds from carbon by high-pressure shock waves. They could make a hexagonal diamond for the first time. In our gray cast iron samples, we also found in graphite nodules hexagonal lonsdaleite crystals with their characteristic Raman triplet at 1244, 1305, and 1356 cm^-1 [1,4]. Sobolev et al. (1993) [1] and Sobolev et al. (2020) [5] synthesized diamonds in gray cast iron; however, they obviously did not check the cast iron for the presence of diamonds before their detonation experiments. So, diamonds/nanodiamonds are probably always present and serve in part as seeds for the formation of larger diamonds. In our case, detonation or shock waves can definitely be excluded.

Samples and Methods

For the study, we selected six different gray cast iron samples, in which we found diamonds. One sample (Sample 2) served for detailed research [6]. Sample 2, a rectangular parallelepiped of 15.2 g, is treated in two steps with hydrochloric acid (25%). In the first step, we solved 1.90 g to remove possible remnants of diamonds from the grinding and polishing. We rejected this part. In the second step, we extracted an additional 1.60 g of iron from the sample. The solid remnants in the hydrochloric acid solution of the second sample (graphite, boron, and carbides) were cleaned with destilled water and placed on a microscope slide to dry. The amount of the remnant was not determined, therefore the exact amount of diamonds is not possible. For the study, we generally used a polarization microscope for transmission and reflection (Olympus BX43) equipped with an X-Y or rotating stage coupled with the EnSpectr Raman spectrometer R532. Details are in Thomas et al. (2025) [7]. Figure 1 shows the reference Raman spectrum of diamond from Brazil (Mining Academy Freiberg: 2453/37) in the range of 0 to 2000 cm^-1, taken at 30 mW on the sample.

Figure 1: Reference Raman spectrum of diamond. The FWHM (Full-Width at Half Maximum) is 4.26 ± 0.42 cm^-1.

Experimental Section and Results

From the second cleaned solution fraction, we put small graphite flakes on the microscope slice. Under many graphite flakes with single or a couple of diamonds, we found two flakes with many diamonds, which are ± regularly distributed (Figure 2). For a more detailed study, we used the number two shown in Figure 2. The diamonds show in graphite a double arrangement of hyperboloids perpendicular to each other, demonstrating two-dimensional stress.

Notes to Figure 2: A picture tells us more than 1000 words. For the distortion of 3.8 µm, a force of 3.8 N is applied over the area of the graphite flake, which corresponds to a pressure of approximately 3.04 GPa. That means, in local areas, high pressure is possible (see further below).

Figure 2: Graphite flake (black) with ± regularly crystallographically arranged diamonds (white). The graphite flakes have a thickness of about ≥ 2.5 µm, and the percentage of diamonds is ~7.4 vol%. The red lines show schematically the bend of the graphite plate perpendicular to the c-axis.

From Raman spectrometric measurements, we obtained the following data for 11 diamonds (Table 1).

Table 1: Results of the measurements on the diamond spheres shown in Figure 2 (532 nm laser, four mW excitation on the sample).

Crystal	Diamond (cm-1)	FWHM (cm-1)	Graphite (cm-1)	FWHM (cm-1)	n
Flake 1b	1329.2 ± 2.8	66.4 ± 0.6	1582.4 ± 7.7	62.2 ± 2.1	11

n: number of studied crystals; FWHM: Full-Width at Half Maximum.

A typical Raman spectrum is shown in Figure 3. Note, according to Prawer (1998) [8], the Raman cross-section for sp₂ clusters is significantly greater (by a factor of 50, including the 1580 cm^-1 band) than that for sp₃-bonded carbon of diamond. The large FWHM values are a problem because all diamonds (in cast iron or in nature: see Thomas, 2025) [9] show such large values for FWHM. It is conceivable that the partial hydrogenation of the diamond surface is its origin [10].

Figure 3: Raman spectrum (4 mW excitation) of diamond from Figure 2. Because the measuring point is about 1 µm in diameter, the Raman band of “graphite” is the result of a relatively high concentration of sp2 cluster in diamond [8], by a thin carbon coating on the diamond, or simply by the excitation of the graphite matrix by chance.

From a rough estimation (Figure 2), we obtain a diamond volume of about 7.4 vol% for the graphite flake. The regular arrangement of the more or less spherical diamond crystals points to a growth of the diamonds in the graphite flake during the production of the gray cast iron. Extreme high pressure to produce diamonds can, in each case, be excluded. The exact origin of such diamonds can not be explained at present. It is clear that the diamonds grow after the formation of graphite, and they are not trapped in graphite by change, as the regular arrangement shows. Obviously, the diamonds mark the grain boundary edges. The low number of diamonds shows that the graphite is less defective and reduces the number of nucleation sites. The process of transforming graphite into a diamond typically requires very high pressures (4.5-6 GPa) and temperatures (900-1,300°C). Such high values for pressure are probably not achieved in the production of gray cast iron. The exact pressure and temperature can vary slightly depending on impurities and the presence of catalysts. Here is iron the favorite. Generally, other catalysts are nickel, cobalt, and alloys of Fe, Ni, Co, and rare-earth metals. These can dramatically lower the pressure and temperature needed for diamond formation. The metals dissolve graphite and facilitate the rearrangement of carbon atoms into the diamond structure at much lower pressure, or even lower with optimized catalysts.

Figure 4 results are approximately obtained using the temperature, the molar volume of graphite and diamond, and the entropy of both using Dean’s (1979) [11] and Jacob’s (1995) [12] data. According to Figure 4, the pressure reduces from 2.19 to 1.82 GPa (triangle), which also corresponds to the stability of calcite as a melt droplet in the gray cast iron (see Brümmer et al. 2025 [6]), which demonstrates a high pressure for the stabilization, because at atmospheric pressure, calcite decomposes into CaO and CO₂ [13]. An increase in the molar volume is, for example, possible by the formation of a porous texture or impurities (Fe₃C). One all-present catalyst is the iron in the surroundings of the graphite flakes, the primary catalyst. Another one is the methane sitting in the graphite lattice, as we assume from Figure 5. This figure shows the sharp Raman forbidden band at 867 cm^-1; however, this band is usually only IR active [14,15]. According to Kawashima and Katagiri (1999) [16], it arises from impurities and/ or structural imperfections. Dmitruk et al. (2012) showed that after quantum-chemical investigations, the methane molecule intercalates between graphite planes. The most stable form of these is methane, which is in the form of dimers or clathrates. The concentration of such intercalated amounts to mole% [14]. With Raman spectroscopy, we can show that most graphite in the gray cast iron contains methane and benzene, which, at high temperatures, behave supercritically and move very fast to the grain boundary edges, forming there clathrates as a preliminary stage of DLC (Diamond Like Carbon).

Figure 4: The graph shows the calculated equilibrium curve of graphite and diamond, as well as the effect of the increase in the molar volume of the graphite (black triangle) by 7.4% resulting in a drop in pressure of 0.95 GPa at 1000°C.

Note: In all studied samples, this remarkable Raman-forbidden line was found, which means that at least this line is typical for gray cast iron. More studies on this problem are necessary. The intercalated methane (CH₄) together with hydrogen (H₂) serves in the supercritical state as seeds for the nucleation of diamonds in the graphite grain boundaries, and with fewer defects.

Figure 5: Raman bands of graphite (D-band at 1350 cm^-1; G-band at 1580 cm^-1). The insertion shows the forbidden Raman line in graphite, which is usually only IR-active. This line is typical for graphite in gray cast iron. The intensity of the forbidden Raman band at 866 cm^-1 is only 1/30 of the G-band.

Discussion

We found diamonds using Raman spectroscopy in all the gray cast iron samples studied. Finding diamonds with classic light- microscopical methods in gray cast iron is somewhat tricky. Raman spectroscopy is very helpful for making a clear decision. To prevent diamonds from being used for the preparation of the samples, we rigorously removed the surface of the samples with hydrochloric acid. During that process, we found occasional flakes of graphite with crystallographically arranged diamonds. Most graphite flakes show only, if at all, a single or a small couple (1-3) of diamonds. In a recent paper by Shumilova et al. (2025) [17], the results of glassy carbon synthesized from supercritical fluid in the C-O-H system at 800°C and pressures of 500 to 1000 atm show similar carbon bodies; however, with strongly different Raman band positions using the 532 nm laser (Table 1 in that paper). Because the FWHM values for diamonds are huge [9], more sophisticated studies are necessary.

Acknowledgement

The author extends his sincere gratitude to Gregor Brümmer (Altzheim, Germany) and Klaus Scheiblauer (Wiener-Neustad, Austria).

References

1. Sobolev VV, Taran Yu N, Gubenko SI (1993) Synthesis of diamond in cast iron. Met Sci Heat Treat 35: 3-9.
2. DeCarli PS, Jamieson JC (1961) Formation of diamond by explosive Science 133: 1821-1823. [crossref]
3. Cowan GR, Dunnington BW, Holtzman A H (1968) Process for synthesizing United States Patent Office, 3,401,019, 12 P.
4. Goryainov SV, Likhacheva Ayu, Ovsyuk NN (2018) Raman scattering in Journal of Experimental and Theoretical Physics 127: 20-24.
5. Sobolev VV, Gubenko SI, Rudakov DV, Kyrychenko OL, Balakin OO (2020) Influence of mechanical and thermal treatments on microstructural transformations in cast ions and properties of synthesized diamond Solid State Physics, Mineral Processing 4: 53-62.
6. Brümmer G, Thomas R, Scheiblauer K (2025) Wie wächst er nun? Der Graphit im Gußeisen unter dem Raman-Spektroskop. Giesserei Rundschau 72, 4: in press.
7. Thomas R, Brümmer G, Scheiblauer K (2025) Unexpected carbon phases in grey cast iron – diamond, calcite, and GEMS 7(4): 1-6.
8. Prawer S, Nugent KW, Jamieson DN (1998) The Raman spectrum of amorphous Diamond and Related Materials 7: 106-110.
9. Thomas R (2025) Diamond in pegmatites and related rocks in the upper GEMS 7: 1-5.
10. Datta J, Ray NR, Sen P, Biswas HS, Vogler EA (2012) Structure of hydrogenated diamond like carbon by micro-Raman spectroscopy. Material Letters 71: 131-133.
11. Dean JA (1979) Lang’s Handbook of 15^th ed. McGraw-Hill, New York. 1291 P.
12. Jacob KT (1995) Determination of the Gibbs energy of diamond using a solid state Solid State Communications 94: 763-765.
13. Hou M, Zang Q, Tao R, Liu H, Kono Y, et al. (2019) Temperature-induced amorphization in CaCO₃ at high pressure and implications for recycled CaCO₃ in subduction Nat Commun 10: 1-8.
14. Dresselhaus MS, Dresselhaus G, Eklund PC, Chung DDL (1977) Lattice vibrations in graphite and intercalation compounds of graphite. Materials Science and Engineering 31: 141-152.
15. Thomas R, Rericha A, Pohl WL, Davidson P (2018) Genetic significance of the 867 cm^-1 out-of-plane Raman mode in graphite associated with V-bearing green Mineralogy and Petrology 112: 613-645.
16. Kawashima Y, Katagiri G (1999) Observation of the out-of-plane mode on the Raman scattering from the graphite edge plane. Physical Review B 59: 62-64.
17. Shumilova TG, Ivanova LA, Isaenko VV, Ulyashev VV, Medvedev Vya, et al. (2025) Glassy carbon synthesis from supercritical fluid in the C-O-H system at 800°C and pressues of 500-1000 Glass and Ceramics 81: 401-408.