The study of regularities in the formation of mineral deposits has always been the main task of geology. By the end of the last century it was found out that the distribution of deposits obeys to certain regularities both in localization in the Earth’s crust and in age. It was realized that changes in the material composition of deposits are due to both the nature of tectonic movements in the evolution of the Earth, and the time of appearance of each shell of the Earth. The greatest success was achieved in the study of tectonic processes in the lithosphere. However, it became clear that it is necessary to clarify the mechanisms and energy of interaction between internal and surface structures. Initially, a one-dimensional model was used to solve this problem. Adiabatic compression and three-dimensional effects due to the deposition of bodies on the surface of a growing planet have not been taken into account. The main source of internal energy of the protoplanet at this stage of its formation was the release of heat accompanying the natural radioactive decay of long-living uranium and thorium. The contribution due to the release of heat from the impact of small particles and bodies of the protoplanetary cloud on the surface of the growing planet turned out to be very small, since most of it, according to Stefan’s law, was radiated into space. As a result, sufficiently low estimates of the temperature of the inner regions of the protoplanet were obtained by the end of the process of their accumulation (Figure 1). As can be seen from these results, the obtained temperature estimates turned out to be lower than the melting point of iron at a given depth. From this it follows that, since the values of the diffusion coefficient of iron in a solid silicate substance are very small, for the entire age of the Earth, the separation of a predominantly iron core from a predominantly silicate mantle could not occur. In order for such separation to be realized, the molten state of the substance and the convective heat and mass transfer of the components in the melt are required. In    , a new model for the accumulation of terrestrial planets was proposed, which uses modern results of isotope geochemical analyzes, which allowed obtaining reliable estimates of the concentration of short-living naturally radioactive isotopes and, above all, 26Al in the matter of the protoplanetary cloud  . On the basis of these data, new estimates of the temperature distribution in the growing pre planetary bodies of the planet in the
Figure 1. The Earth’s temperature distribution after the end of accumulation. 1―the temperature of full mantle material melting; 2―the temperature of the solid mantle; 3―the melting iron temperature; 4―the temperature distribution into the Earth by Safronov  (lower curve); 4―the temperature distribution into the Earth by Lubimova (upper curve)  .
“feeding” zone of the Earth were obtained. It was shown in  that the energy release in the decay of short-living radioactive elements in small bodies is sufficient for the temperature inside such a protoplanetary body to become greater than the melting temperature of iron. This ensures the realization of the process of substance differentiation in density and the development of convection in the inner shells. Originally, a one-dimensional mathematical model was used to numerically solve problem (1)-(4), based on new data on the content of 26Al in the pre-planetary matter of the Earth.
The interest of geologists in the problems of geodynamics, the development of the structure of the interior regions of the planet has increased significantly after the publication of the paper  . The novelty of the presented work is the construction of a quantitative description of the earliest stage of the evolution of the Earth, starting with the process of accumulation from the protoplanetary cloud and tracing the development and development of jet structures (plumes) in the mantle and their manifestation on the surface. This became possible only on the basis of a numerical solution of the problem for a system of nonlinear equations in a three-dimensional model with thermal and real heterogeneities and taking into account phase transitions    . The results obtained, in contrast to   , made it possible to describe the early formation mainly the iron core of the Earth for the first 30 - 50 million years. In comparison with the very advanced results of  , we did not encounter here the necessity, characteristic of “N-body problems”, of setting unknown boundary conditions on the boundaries of the chosen domain of the solution of the problem. In contrast to   , in the realization of the solution of the problem in the three-dimensional model, it was possible to take into account, on each time layer, a random distribution of the thermal and real heterogeneities created by the precipitation of bodies and particles upon impacts on the growing surface of the Earth.
2. Statement of the Problem and the Mathematical Model
The solution of the problem should describe the formation of a protoplanet from the matter of bodies and particles of a protoplanetary solar cloud. The statement of the problem is described in detail in     . The process of accumulation of the planet and its interior regions was accompanied by the formation of primary mantle thermal and material heterogeneities. Thermal heterogeneities are due to the difference in the kinetic energy of the collision of the accumulated bodies with the growing protoplanet, and the real heterogeneities are caused by the difference in the composition of the drop-out bodies and their melting points. To obtain a mathematical description of the process, a boundary value problem is formulated for a system of equations describing the growth of the protoplanetary mass (the Safronov equation is used); balance of momentum (the Navier-Stokes equation); Stefan’s problem for describing the motion of the boundaries of the regions of phase transitions. On the surface of a growing planet, a condition is given, taking into account the random distribution of bodies over the surface and the kinetic collision energy. The 3D model is used. In    , a model for the accumulation of terrestrial planets was proposed, which uses modern results of isotope geochemical analyzes, which allowed obtaining reliable estimates of the concentration of short-living naturally radioactive isotopes and, above all, 26Al in the matter of the protoplanetary cloud  . On the basis of these data, estimates of the temperature distribution in the growing pre planetary body of the planet in the “feeding” zone of the Earth had been obtained in   . The purpose of this part of the work is to obtain a three-dimensional temperature distribution and trace the development of thermal heterogeneities in the mantle at the stage of formation of the Earth, which subsequently provided heterogeneities of mineral deposits localization. To calculate the growth velocity of the mass of the arising planet, Safronov model describing the mechanism of planet accumulation is used  :
where: ω is the angular velocity of the orbital motion, σ is the surface density of matter in the “feeding” zone of the planet, m(t) is the pre planetary mass, M is the modern mass of the planet, R is the radius of the growing body, θ is the statistical parameter taking into account the particle distribution by masses and velocities in the “feeding” zone, m0, R0 is the mass and radius of the arising planet at the initial moment of time. For each achieved mass value, the distribution along the radius of density, hydrostatic pressure, and melting point in the inner regions is calculated. The density distribution in the mantle is calculated from the Williamson-Adams equation, which can be written as:
where ―density, ―density increment, ―radius increment, ―
density on the Earth’s surface, ―gravitational acceleration. Here G
is the gravitational constant, m is the mass enclosed inside the sphere of radius r. ―a seismic parameter that can be defined as:
where K is the modulus of compression. The distribution of the hydrostatic pressure along the radius is found from the Euler equation describing the mechanical equilibrium of a liquid  .
where p(r) is the pressure. In the core, mainly of iron composition, the dependence of the melting point on pressure is calculated from  . In the mantle formed mainly by silicates, the dependence of the melting point on pressure is used, as suggested in  . The temperature at the surface of the layer with the power Δr formed during the accumulation of the planet during the time Δt is calculated from the equation ensuring the balance of the incoming part of the potential energy of the gravitational interaction of the bodies, the heat input for heating the incoming material and re-emitted into the heat flux space, taking into account the transparency of the external medium   :
where: is the density of matter, G is the gravitational constant, m is the mass of the growing planet, R is its radius, T is the temperature at the outer boundary, T1 is the temperature at the sunflower point, ε is the transparency coefficient of the medium, cp is the specific heat, k―the fraction of the potential energy converted into heat. The temperature distribution in the body of the increasing radius is found from the solution of the heat conduction equation with allowance for the possibility of a melt appearance without explicitly isolating the position of the phase transition boundaries, convective heat transfer in the melt   :
Here is the temperature at time t, ―here is the mass velocity of the substance, ―the capacity of the internal heat sources, the effective heat capacity and thermal conductivity, which take into account the heat of fusion in the Stefan problem  and the presence of convective heat transfer, TR is the temperature on the surface of the growing planet, obtained from Equation (5), T0 is the initial temperature distribution, is the Nabla differential operator. The value of the effective heat capacity was calculated in accordance with the idea of the through-account method  , a delta-shaped function for the specific heat was constructed such that in the temperature interval ΔТ at the phase transition point the effective heat capacity was:
where L is the heat of the phase transition. Convective heat transfer in the melt layer is accounted for by the effective coefficient of thermal conductivity   :
where Ra is the Rayleigh number, Rak is the critical value of the Rayleigh number; α is the coefficient of thermal expansion, h is the thickness of the convective layer; ΔT is the temperature difference at its boundaries; g is the magnitude of the gravitational acceleration; v is the effective viscosity in the layer. The regions of the melt were compared with the melting point and temperature, for a given moment of time. The power of internal heat sources is the sum of the heat release during radioactive decay of 26Al and adiabatic compression:
where is the concentration at the initial time; λ is the decay constant; q is the energy of decay; 26Al/27Al is the isotope ratio  , Q―summarized capacity of the inner thermal sources, Q1―capacity of heat release by adiabatic compression  , NA―number of Avogadro, ―concentration of Al2O3, ρ―density.
The boundary value problem was solved numerically. When solving Equation (5), using a random function that is a standard random number generator, we get the distribution of thermal heterogeneities in terms of kinetic energy and over the surface of the growing planet. The boundary value problem for Equation (6) written in spherical coordinates was solved by the method of finite differences using the splitting scheme   . Steps on the spatial and temporal grid are no uniform.
It is thermal and material heterogeneities that are the objects, which are registered by the methods of modern geophysics and geochemistry in the mantle. As can be seen from the presented results, they are formed during the accumulation of the Earth and form heterogeneities of various sizes and shapes. The resulting inhomogeneous inclusions change with time. In the presented model variants of their evolution are considered. which are obtained as a result of a numerical solution of the boundary value problem for a system of differential equations describing the change in the mass of a growing planet using the Safronov equation, the flow of matter using the Navier-Stokes equation, the motion of phase boundaries using the Stefan problem; On the surface of a growing planet, conditions are formulated that take into account the random distribution of bodies and particles as a function of the energy and the point of incidence released. Since this distribution of falling bodies is given randomly, and it randomly changes during the accumulation of the model of a growing planet, it was received a set of different solutions to the problem. As it shown by the different results of the solution, to the end of the planet’s accumulation a relatively thin heterogeneous crust forms, and in the upper mantle an asthenosphere layer. In the mantle, heterogeneous inclusions are observed, ranging in size from the first to several hundred kilometers.
The insets show the results on an enlarged scale. Radial chains of heterogeneities of the elevated temperature in the mantle are seen for the instantaneous time.
As it is shown in the upper inset to Figure 2, regional molten temperature heterogeneities can be traced to depths of about 1 km from the surface of the planet at the final stage of its formation. A more detailed further behavior of the molten regions in the variants obtained could not be isolated due to the spatial resolution of the realized numerical grid. The obtained results make possible to trace the change in the thermal structure of the outer core and mantle as accumulation and growth of the Earth’s radius. In connection with the considered problem of the evolution of the processes of metallogeny, special attention should be paid to the traceable unification of local thermal and material heterogeneities and the active formation of increasing structures taking on an increasingly vertical form. In the mantle are observed heterogeneous inclusions, ranging from the first to several hundred kilometers with a complex interface. These
Figure 2. A variant of the results of numerical modeling of the dynamics of thermal heterogeneities toward the end of Earth’s accumulation.
inclusions under the action of the Stokes force begin to float, leading to the formation of jets or so-called mantle plumes. Floating local thermal inclusions and formed jets reduce the convective stability of the mantle and promote the development of regular convective currents. The very dynamics of heterogeneous inclusions, their random distribution in temperature and material composition, provided a special metallogeny of the cratons, which no longer repeated in the geological history of the planet. In particular, it caused the fact that contrary to the expected, the magmatism belonging to the effect of plumes of the same age and depth of embedding is not sufficient to form a close material composition arriving to the earth’s surface and this can be a consequence only of contamination processes.
On Figure 3, in contrast to Figure 2, the results obtained for the same models as those shown in Figure 2, are presented in the 3-dimensional segment. The distribution of 3-dimensional temperature inclusions is shown for the time moment close to the completion of the planet’s growth. It can be seen that by this time long jets or plumes extended along the radius are formed.
The numerical solution of the formulated problem for a three-dimensional model with allowance for the random distribution of bodies and particles to the surface of a growing Earth for the first time made it possible to trace the further fate of thermal heterogeneities in the nucleus, mantle and crust in formation. The thermal conductivity of a large number of primary, due to the composition of the trapped bodies of the protoplanetary cloud, increased in comparison with
Figure 3. Distribution of temperature and thermal heterogeneities to the completion of Earth’s accumulation in the 3-dimensional sector. Light-blue chains of heterogeneities, more heated than the surrounding matter of the mantle, trace the formation of primary plumes, some of which then rise into the emerging crust.
the host material, leads to the formation of jets, plumes, and the differentiation of accumulated matter in them in the gravitational field. Unlike the further formation of jets in the modern Earth structure, when they are formed mainly on the core-mantle section and the 680 km phase boundary in the mantle, these plumes could be formed at arbitrary depths depending on the size and position of the primary heterogeneities. Thus, at this stage of the development of the Earth P-T, the conditions for separation of the magmatic melt and its composition could be extremely diverse; this is the reason for its difference in composition from the modern one. This requires a careful study of the dynamics of the emerging binary system Earth-Moon. Since in the process of its formation, the distribution of pressure and temperature in the inner regions of the emerging Earth, the velocity of the planet’s orbital rotation may substantially differ from these thermodynamic parameters in the Earth as alone body.
These results had been achieved by support of the grant RFBR № 16-05-00540.