Kaula’s rule – a consequence of probability laws by A. N. Kolmogorov and his school

RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 19, ES6007, doi:10.2205/2019ES000651, 2019

Kaula&s rule — a consequence of probability laws by A. N. Kolmogorov and his school

E. B. Gledzer1 and G. S. Golitsyn1

Received 18 January 2019; accepted 28 November 2019; published 12 December 2019.

The paper analyzes a possible cause for the universal behavior of the covariating fluctuations of planet gravity. The consideration based on the idea that the topography fluctuations are governed by a random Markov process leads to universal dependence fc-2with k being the amplitudes of the j-th spherical harmonic of the terrain profile. This law known as the Kaula&s rule is then derived from the solution of the Fokker -Plank equation for the fluctuations of the terrain profile as the function of the horizontal coordinate. The respective diffusivities for Earth and Venus are retrieved from the existing data. KEYWORDS: Kaula&s rule; Thumb rule; planet gravity; Fokker-Planck equation.

Citation: Gledzer, E. B. and G. S. Golitsyn (2019), Kaula&s rule - a consequence of probability laws by A. N. Kolmogorov and his school, Russ. J. Earth. Sci., 19, ES6007, doi:10.2205/2019ES000651.

Introduction

The Kaula&s rule [Kaula, 1966] sometimes called "the rule of thumb" exists for over half a century. However, its origin is unknown. In its original form it claims that the terrestrial gravity field fluctuations expanded in spherical harmonics have the spectral asymptotics proportional to j-2, where j is the number of the harmonic. Similar behavior was established to hold on Moon, Mars [Rexer and Hirt, 2015], and Venus [Turcotte, 1997]. In this decade the same property has been found for asteroid Vesta whose diameter does not exceed 300 km [Konopliv et al., 2014] and for asteroids of few kilometers in size [McMahon et al., 2016]. In later 1980-s it was found that surface reliefs of Earth and Venus have the same j-2 asymptotics [Turcotte, 1997]. In [Rexer and Hirt, 2015] the harmonics with j up to tens thousands have been computed. This means that gravity field fluctuations and relief are connected. Therefore "the rule of thumb"

Copyright 2019 by the Geophysical Center RAS. http://rjes.wdcb.ru/doi/2019ES000651-res.html

is expected to be a law whose nature still awaiting for explanation.

Here we use the Fokker - Planck equation in the form presented by A. N. Kolmogorov, [1934] and its further results of its statistical analysis developed by his pupils A. M. Obukhov [Obukhov, 1959] and by A. S. Monin and A. M. Yaglom, 1975]. The basic equation formulated in [Kolmogorov, 1934] is more general than it is necessary here. Later this equation was used for studying other natural phenomena: like turbulence in [Obukhov, 1959], and recently sea wind waves spectra, the cumulative distribution of litospheric plates, etc. in [Golit-syn, 2018]. This impelled us to return to an attempt to derive analytically the Kaula&s rule. For reader&s convenience we present the description of ideas and methods of [Kolmogorov, 1934; Obukhov, 1959; Monin and Yaglom, 1975] rarely used in general geophysics but still appearing in various papers and books during past decades.

Basic Equations

In the beginning of 1930-s A. N. Kolmogorov had published three papers on analytical methods in

the probability theory. The two-page work [Kol-mogorov, 1934] contains the essence of his approach started by A. Einstein and developed further by Fokker and Planck. Kolmogorov proposed a fundamental solution to the evolution equation for the probability distribution p(ui, Xi, t) of the 6D vector u, x provided a delta correlated force causes the respective Markov process. This is an approximation of the processes where correlation time of the random forces is much shorter than the reaction time of the considered system. In 1D approximation the action of respective random accelerations at with { a(ti)a(^2)) = a2S(t i — £2) can be described by the simplest equations [Gledzer and Golitsyn, 2010],

dt dt where u( ) and X( ) are the velocity and displacement of a particle under random acceleration a(t).

A. N. Kolmogorov used a Fokker-Planck type equation [Kolmogorov, 1934] to describe the evolution of the probability density p(ui,Xi, t),

dt ° dx

dp ^ d2p du2 &

where D is a coefficient of diffusion in the velocity space. The fundamental solution2 to eq. (2) has the form (see also [Monin and Yaglom, 1975], §24.4):

p(ut)={2Ê# ) (

D t2 & D ts

A. M. Obukhov [Obukhov, 1959] was the first who analyzed and used this equation. He showed that the coefficient D is proportional to the generation/dissipation rate of the kinetic energy of the random motions D = e/2. The details of his analysis are presented in [Monin and Yaglom, 1975; Golitsyn, 2018]. The solution (3) demonstrates that the seeked probability distribution has the normal character. This solution has two scales (angular brackets mean ensemble averaging):

£ t3,s = 2D.

Figure 1. Temporal behaviour of the second velocity and coordinate moments for different particle ensemble sizes as in [Gledzer and Golitsyn, 2010] (see text).

Combining the scales for velocity and coordinates Obukhov formulated the scale for the mixed moment { UjXj) = K in terms of £ as

K = £t2. (4&)

The time dependencies eqs (4) and (4&) were verified numerically for the ensembles of N randomly moving particles [Gledzer and Golitsyn, 2010]. Figure 1 demonstrates that even for N = 10 the relationships (4) fulfil satisfactorily as it is shown in Figure 1. Moreover, on introducing the dimen-sionless variables to non-dimensional ones (denoted with tilde overhead)

t = tr,u = U(sr)1/2,r = r(srà)

reduces eq. (2) to the universal or selfsimilar form [Gledzer and Golitsyn, 2010]

dp dp —- +1-—

The scales (4) can be interpreted as the structure functions with zero initial conditions which allow us to introduce their spectral forms [Monin and Yaglom, 1975]. If we replace t by r (see eq. (4)), exclude time from the second scale (4) and put it into the first scales and to (4&) and substitute the result into eq. (4&), we obtain the results of Kolmogorov - Obukhov of 1941 Du(r) = (er)2/3 and Richardson - Obukhov [Obukhov, 1941; 1959;

Richardson, 1926; 1929] for the eddy diffusion coefficient:

= Л/3_4/3

These are the results for the inertial interval of turbulence. The results as relations (4) are valid irrespective of the space dimension and the numerics prove it [Gledzer and Golitsyn, 2010]. It is surprising that since 1920s nobody tried to understand why the Richardson&s "law of- 4/3" holds up to 2-3 thousands km in the horizontal direction [Richardson, 1926; 1929], i. e. the vertical component and full velocity isotropy are not important. In 1999 Lindborg [Lindborg, 2008] analyzed the experimental data of commercial flights in terms of structure functions. He showed that the dependence r2/3 is traced up to distances of 2 thousands km. His results in terms of the eddy diffusion explain the Richardson results of 1920-s. These facts can be explained starting with a two-page note of Kol-mogorov, 1934 [Kolmogorov, 1934], and the following development of his results by Obukhov, [1959] and Monin and Yaglom, [1975]. Only the Markov properties of the process are important rather than the full isotropy except in horizontal direction. Our numerics of 2010 [Gledzer and Golitsyn, 2010] confirmed this independence of the space dimensionality for the second moments with the accuracy of the numerical coefficients up to 0 (1).

Spectral Representations

For further development we need an analytical connection between the second moments and their spectral representations when both of them are power laws. Let us consider first the structure function D(t) = At7, 0 < 7 < 2. The relation to the spectrum is [Monin and Yaglom, 1975]

D(t) = 2 (1 - cos wt)E(w)dw,

,C = -Г(7 + 1) sin 7, (8) ■к 2

where r is the gamma function. The length scale related to the mean square displacement of a randomly moving particle has therefore the spectrum

Сш-4, С = А/ж and for the velocity scale (u2) = et the spectrum is eu-2. According to the terminology introduced by Yaglom, [1955; 1986] the spectrum ш-2 corresponds to a process with stationary increments of the first order (0 < ^ < 2), while the spectrums ш-4 has the stationary increments of the second order (2 < ^ < 4). In order to avoid the singularity at ш ^ 0 the transformation kernel in eq. (7) should be (1 — coswt)n, with n = 1 for velocities, and n = 2 for coordinates.

Toward Thumb Rule

When a flying altimeter measures the relief, we obtain a temporal variations of height signal h(t). The time t here is related to horizontal coordinate as y = ut, u being the flight velocity. The standard Fokker - Planck equation with stochastic velocity fluctuations in this case is

dp p д 2 p dt = dh2 If y = ut it transforms to

dp = D d^p D = D dy dh2 & и&

The diffusion coefficient D\\ has the dimension of length. But if the vertical Lz and horizontal Ly scales are different, then ^ = L2,/Ly which reflexes clearly the diffusion character of the relief formation process during millions years. Equation (9) yields the scale and the structure functions for zero initial condition in the form:

(h2(y)) = Diy with the corresponding spectrum

S (к) = Dik-2, к = 2ж/Х

The inversion gives for the spectrum [Monin and Yaglom, 1975; Yaglom, 1955; 1986]:

Xy being the corresponding horizontal wave length. This equations applicable for small scale continuous areas. The spectrum of the relief slopes in this approximation is S^(k) = k2S(k) = D = const. The constancy of a spectrum is referred to as "white noise". The slope of the relief

dh(y) dy

= Ç (У)

for which we accept the hypothesis of ¿-correlated horizontal coordinates

( ç(î/i)ç(î/2)) = (y i - V2)

= ft ( V, 0)

is the slope angle. Let us introduce x = cos 9 and present the relief as a sum of normalized associated Legendre polinomials P1™^

Equation (12) corresponds to the Markov character of forces in [Kaula, 1966; Rexer and Hirt, 2015; Turcotte, 1997]. Here s2 may be interpreted as the mean square slope angle, a dimensionless value.

Now we return to the reliefs of concrete planets. Recently an extensive paper on ultrahigh spherical harmonics analysis for topography of Earth, Mars and Moon has appeared [Rexer and Hirt, 2015] with the horizontal resolution up to hundreds of meters. Much earlier analyses of topography for Earth and Venus were given in [Turcotte, 1997]. In that book there are also the data for mountain, hilly and plain areas of hundred kilometers in size at the Oregon state. The flights were done in various directions in each type of terrain. For records in the range 1-60 km the mean exponent value of these spectra were found [Golitsyn, 2003] to be 2.03 ± 0.04. Unfortunately, there were no data in [Turcotte, 1997] for spectral amplitudes. Obviously k-2 small size spectrum is a continuation of what Kaula found for large scale harmonics.

In the book [Turcotte, 1997] the spectral exponent 2 (and other statistical measures, such as Hausdorff measure, the fractal index d) were retrieved from the original records. But no explanation of particular values was proposed.

The analysis done by Kolmogorov [Kolmogorov, 1934; Monin and Yaglom, 1975] several decades earlier was clearer and simpler. The only assumption on -correlation of slopes applies for the processes in consideration, i.e. the Markov&s hypothesis. Such assumption is used in many branches of theoretical physics meaning only that the system reaction time is much larger than the correlation time of stochastic influence on it [Golitsyn, 2003].

Spherical Harmonics

For large planetary scale features the spherical harmonic analysis is used. Let us consider the relief on a sphere z( f, d) with a radius r, where f is longitude, 0 < f < 2k, d is a co-latitude, 0 < d < k. Then

dz(f, 9)

z( p, d)

= -exp(^mvHml(x), (13)

j=1 m=-j

(x) = N^P^ (x),

f 2j + 1(j -M)! V/2

V 2 (j + M )!J &

Then due to eq. (13) the relief slope is expressed through the derivative,

. d$&H (x)

(e(0, 0) = E E a™exp(2mip)sindj=1 m=-j

Consider the mean normalized energy (in quadratic sense) of the relief slopes as follows,

f27T f1 2 Ee = J J < (22 > dpdx =

j=1 m=-j

f\\ . (x) \\

j-i(l -x2)(-*-)dx

We consider the mean values of the spectral components < |am|2 > to be independent of m, which is confirmed by observations and on practice [Kaula, 1966; Rexer and Hirt, 2015; Turcotte, 1997]. It thus depends only meridional index j and can be found on summing over longitudinal index m:

< Kj2 >=

¿0&+ 1)0&+ 1/2)& where a is the mean angle of the relief slope. Taking this into account and calculating the integrals in eq. (15) we have on summing over m

Ee = 2^ESi, 3=1

We thus obtain the constant spectrum of "white-noise", i.e. -correlation of Then for the relief energy (as mentioned above), we derive a quadratic value integrated over the spectrum from (13),

E = J J < z2 > dipdx = (18)

œ J /»1

j-1 m=-j

E = 2ir£ E

j=1 m=-j

r2a2 (2j + 1)

=1 i0& + 1)0& + 1/2)

Here D-t is the horizontal diffusion coefficient and

E KJ+I) = 1 <"ote 5(5+1)

idenp\\™\\(cos 0) « (-W J \\ij cosO

/ - mi \\ f. 1 \\ _

x co^a - 4 + ^J , « +

Now instead of co-latitude we will use the length along the meridian assuming = 0 at the north pole,

y = r0 < d <i.

Then we rewrite the phase a in eq. (20) as a = ydk(j + 1), where dk = ^ is an increment of the wave number k = j ■ dk at j ^ x. Then for very small horizontal scales one has at dk ^ x,

a = ky, k = j ■ dk = j/r

j-1 m=-j

It is clear from here that the mean of the spectral components (17) provides the spectrum relief ~ j-2 at j >> 1, as was found from the empirical data of refs [Kaula, 1966; Konopliv et al., 2014; McMahon et al., 2016; Rexer and Hirt, 2015; Turcotte, 1997].

Now let us apply the wave number k in the eq. (18) for the relief energy. At j ^ x we obtain

E = 2iD1k-2 = 4ira2k-2.

tity). The dependence of our harmonic coefficients 1/i(i + 1) (see eq. (19)) differs from Kaula&s 1/j2 by the factor 1/j2(j + 1). At j = 4 it is equal to 1/80.

In order to replace integers by continuous wave numbers j at j ^ œ, i.e. for small horizontal scales, we shall use in eq. (13) the asymptotics of the polynoms p]m& (cos d) at j ^ œ :

Here 4^ra2fc-2 = S(k) is the spectrum with dimension of the cubic length (recall that a the rms slope angle). It can be seen that the value ra2 in eq. (23) is proportional to the diffusion coefficient D1 in the coordinate space.

The relief spectrum in eq. (19) differs from the Kaula&s "rule of the thumb" a j-2 by the multiplier [j(j + 1)]-1. The relative difference between the two at j ^ x decreases as [j2(j +1)]-1 but for small j it is quite noticeable, e.g. at j = 1 the difference is 1/2, and at j = 4 it is 20%. This shows that the hypothesis on the "white noise" of the relief slope angles is not accurate at these scales. Thus the corresponding relief structure at j < 5 could be different for each celestial body due to its peculiar tectonics.

Results and Discussion

Our scheme of the relief evolution basing on eq. (9) of the probability density changes is the horizontal diffusion of the tectonically formed vertical structure under the gravity acting along slopes. Slopes as a measure of forces resist to winds; water and rocks are running down along them, etc. Figure 2 (Figure 7.19 from [Turcotte, 1997]) presents the amplitudes of the relief spherical harmonics: for Earth j < 180 and for Venus j < 50. At j > 4 all coefficients approach to the inverse square dependence.

We compare our theoretical expression for harmonic coefficients from eq. (19)

Figure 2. Relief spectral harmonics for Earth and Venus after [Turcotte, 1997].

related only to equatorial plates [Turcotte, 1997]. The global relief characteristics can be estimated with much higher precision and efforts (see e.g. [Rexer and Hirt, 2015]), which we leave for much younger colleagues, but we believe that it would be worth doing.

Figure 3 presents the relief spectra as the functions of the surface linear wave length k. The authors of [Rexer and Hirt, 2015] have computed spherical harmonics for Earth and Moon up to = 46200 and for Mars up to j = 23100. Our eqs. (20) - (23) present high frequency harmonics in terms of k-2. At Figure 3 we added straight lines corresponding to the spectra k-2 and k-4, the last for very high wave numbers. The dependence k-2 holds well for Mars, line 1, and for three terrestrial spectra, lines 3-5, with some steeping for waves shorter than a few kilometers. The Moon spectrum falls down and remains almost flat at 200 - 80 km, being very specific in this respect. This evidences that the Moon is bombarded by the asteroids the size of about 10 km, since the size of a crater is 10-20 times exceeds the diameter of the asteroid.

The explanation for the relief spectral decrease from k-2 to k-4 can be explained by rejecting the -correlation for the relief slopes hypothesis. The simplest correlation function is

Be(%i)0(V2)) - exp(-%), y = m - yi, (25)

where 3 1/h and h is mean relief height. Its spectrum is

Se (k) =

k(/32 + k2)

= 4irrDi

S = KTTI)

and if we know the spectrum of the derivatives, i. e. angles, the spectrum of the value itself, relief, would be

to their empirical values from Figure 2.

The comparison has been done for selected harmonics with accuracy about 10% within their value order. The harmonic numbers for Earth were 5, 10, 20, 30, 40, 60, 90. This produced Di = 1.3 ± 0.3 m. For Venus the numbers 5, 10, 15, 20, 30, 40 and 50 produced D1 = 0.19 ± 0.03 m. Our diffusion coefficient is determined with accuracy about 20%. The striking difference between these two planets has been already noted by [Turcotte, 1997]. However, one should recall that Venusian data are

Se (k) = k2S2e (k) =

n(32 + k2)&

This spectrum has both right asymptotes: at h2 ^ 32 it would be proportional to k-2 and in inverse case for small distances it — k-4. The gravity acceleration on the Moon is six times weaken that on Earth and other conditions equal the relief inhomogeneities there could be formed easier and higher.

Figure 3. Relief spectra as functions of the surface wavelength for Moon, Mars and Earth [Rexer and Hirt, 2015](3 cases labeled at lower left). We added lines fc-2 and k-4 for large and small distances.

Conclusion

Our analysis is based on the Fokker - Planck -Kolmogorov, FPK, equation cast into eq. (9) for altimeter records. It has the solution (h2) = D1y for the mean square of the relief height. Here D1 is the horizontal diffusivity of the relief. Such a solution has the spectrum S(fc) = Dk-2, with fc = 2i/X along the horizontal direction. The modified FPK equation is based on the hypothesis of the time S-correlation of acting forces which in this case corresponds to horizontal ¿-correlation of the relief slope angles. The spherical harmonics of the relief have been found by Kaula to follow the inverse-square rule for the harmonic number j. Our "white noise" hypothesis on the horizontal spectrum of the slopes slightly modifies the Kaula&s rule for high harmonics by replacing the inverse square dependence by the factor [j(j + 1)]-1. The difference decreases as [j2(j + 1)]-1. The numerator in eq. (24) gives the constant for the Kaula&s rule.

Rather crude estimate of diffusion coefficients D1 performed for Earth yields D = 1.3 m, several times higher than for Venus. A similar investigation may also be performed for the planetary gravity field fluctuations with similar results, though in this case one should account for the internal

structure and tectonics. Consideration of statistical scales (4) as structure functions with zero initial conditions allows us to propose their spectral representation. For the very high harmonics with continuous spectrum the spectral behaviour changes from fc-2 to k-4, see eqs. (25) - (27).

It can be concluded that the Kaula&s rule for gravity field as well as for the relief is an asymptotic consequence of the random walk laws, for the first time formulated by Einstein (see [Bridgman, 1921]) and in the most general and practical way by Kolmogorov and his successors in 1934 [Monin and Yaglom, 1975; Obukhov, 1959; Yaglom, 1986].

Thus the line of study may be considered as the model, the base for use the dimensional analysis for investigation of various processes. When obtaining dimensional formulas one always should remember that the formulas must contain a numerical coefficient which could be retrieved from the data of observations. The formulas can be valid only in a certain range of dimensionless similar numbers. In certain situations there can be more than one such number. How should one act in these cases? One should consult to professor G. I. Barenblatt, 1996; 2003], who was Kolmogorov&s student.

Geophysical interpretation of these results is waiting for a thorough analysis. Because the Kaula&s constant for gravity field is used to obtain the information on the internal structure of celestial bodies (mascones, isostasy), the diffusion coefficient D1 could contain the information on tectonics, surface material properties, etc. This coefficient is an analogue for the Kaula&s gravity field fluctuation constant. Much work has yet to be done by much younger people since the total age of the two present authors is not far from 160 yrs.

Acknowledgment. We thank A. A. Lushnikov for carefully reading the manuscript and indicating some errors in the formulas. This work has been particularly supported by the RAS Presidium Program 7 "The development of non-linear methods in theoretical and mathematical physics". The authors are very thankful to Dr. O. G. Chkhetiany for discussion of the results and his continuous help during the long work.

References

Barenblatt, G. I. (1996), Scaling, Similarity and

Intermediate Asymptotics, 386 pp. Cambridge Univ. Press, Cambridge, UK. Crossref Barenblatt, G. I. (2003), Scaling, 171 pp. Cambridge Univ. Press, Cambridge, UK. Crossref Bridgman, P. W. (1932), Dimensional Analysis,

Kaula, W. M. (1966), Theory of Satellite Geodesy,

Ann. Math. , 35, 166-117. Crossref Konopliv, A. S., et al. (2014), The Vesta gravity

field, Icarus, 240, 103-117. Crossref Lindborg, E. (2008), Stratified turbulence: a possible interpretation of some geophysical turbulence measurements, J. Atmos. Sci., 65, 2416-2424. Cross-ref

McMahon, J. W., et al. (2016), Understanding the Kaula&s rule for small bodies, Proc. 47th Lunar and Planetary Sci. Conference, AGU, Washington, DC. (https://www.hou.usra.edu/meetings/lpsc2016/pdf 2129.pdf)

Monin, A. S., A. M. Yaglom (1975), The Statistical Hydromechanics, vol. 2, §24.4, MIT Press, Cambridge. (Russ. ed. 1967) Obukhov, A. M. (1941), On the energy distribution in the spectrum of turbulent flow, Proc. (Doklady) USSR Ac. Sci., 32, No. 1, 22-24. Obukhov, A. M. (1959), Description of turbulence in terms of Lagrangian variables, Adv. Geophys., 6, 113-115. Crossref Rexer, M., C. Hirt (2015), Ultra-high spherical harmonics analysis and application to planetary topography of Earth, Mars and Moon, Surv. Geophys., 36, No. 6, 803-830, Crossref Richardson, L. F. (1926), Atmospheric diffusion shown on a distance-neighbour graph, Proc. Roy. Soc., A110, No. 756, 709-737. Crossref Richardson, L. F. (1929), A search for the law of atmospheric diffusion, Beitr. Phys. für Atmosphäre, 15, No. 1, 24-29. Turcotte, D. L. (1997), Chaos and Fractals in Geology and Geophysics, 2-nd ed. , 416 pp. Cambridge University Press, New York. Crossref Yaglom, A. M. (1955), Correlation theory for processes with stochastic stationary increments of the nth order, Math. Collection, 37, No. 1, 259-288. (in Russian)

Yaglom, A. M. (1986), Correlation Theory of

Stationary and Related Random Functions, 1. Basic Results, 526 pp. Springer Verlag, Berlin.

Corresponding author:

G. S. Golitsyn, A. M. Obukhov Institute of Atmospheric Physics, RAS, Moscow, 119017, Pyzhevsky 3, Russia. (gsg@ifaran.ru)