A Method for Quantitative Interpretation of Stationary Thermal Fields for Layered Media

A new method to solve thermal conjugacy problems is presented for layered models with a thermal conductivity jump at their boundaries. The purpose of this method is to approximate the inverse thermal conductivity coefficient, which has breaks, by using a combination of step functions. A generalized continuous operator is constructed in a continuous space of piecewise–homogeneous media. We obtained an analytical solution for the stationary problem of heat conjugacy in the layered model with finite thickness and with Dirichlet–Neumann conditions at the external boundaries. An algorithm was constructed for downward continuation of the heat flux to depths that correspond to the top of the mantle layer. The advantages of this method are illustrated by testing the crustal seismic, gravity and geothermal data of a study area in the Urals and neighboring regions of Russia. We examined statistical relations between density and thermal parameters and determined heat flux components for the crust and the mantle. The method enables a downward continuation of the heat flux to the base of the upper mantle and allows us to determine the thermal effects of the lateral and vertical features of deep tectonic structures.


Introduction
Heat flux is an important character of the thermal state of the lithosphere. It indicates conductive heat exchange in solid bodies and defines the direction of the evolution of deep tectonic processes. The use of geothermal data as a part of a complex interpretation scheme (along with gravity, magnetic or seismic data) lowers the ambiguity of geologic modeling and increases the reliability of the obtained models [1][2][3].
Heat flux is measured by temperature and its gradients at deep and ultra-deep parameter wells. These measurements have unique features. First of all, it is necessary to drill in order to obtain direct measurements, which significantly increases survey time and resources. In addition, a thermal equilibrium should be attained in the borehole, and the duration of this process may be up to several months [4]. The data obtained from the well are averaged, and then the analytic continuation of these heat flux values to the Earth's surface is performed for all of the wells. The resulting values constitute the observed heat flux.
Obtaining a measured data set is very difficult because it is only possible to perform surveys of thermally stabilized wells. Thus, any three-dimensional (3D) model is an interpolation of several one-dimensional (1D) data sets, which are measured in such wells. This is why many researchers tend to utilize data obtained by other geophysical methods (such as seismic and gravity data) for thermal model construction [2,3]. One-dimensional measurements do not allow for the detailed study of the volume Poisson integral type and can be calculated numerically. This is one of the reasons why gravithermal models are widespread in theoretical and practical geophysical research [5].
In the case of non-homogeneous thermal conductivity, a refracted component of the thermal flow appears, and the heat is transferred along the path of least thermal resistance. In piecewise-homogeneous models, abrupt changes in thermal conductivity are allowed. The normal component of the temperature gradient at a thermal contact boundary also undergoes abrupt changes. The contact boundary of media with different thermal conductivity values is the place of material discontinuity. A break in normal temperature gradients forms the equivalent simple layer of inducted sources [15,16].
A determination of heat transfers across this simple layer is obtained from boundary conditions of the thermal conjugacy (conditions of the fourth order). The problems of thermal conjugacy may only be solved by numerical methods using a conservative difference scheme for keeping the energy balance for every discrete element [14]. This problem is algorithmically difficult and requires a significant amount of resources, even for two-dimensional cases [8,9].
We developed a new analytical method for the thermal conjugacy formulation as a boundary value problem for the potential fields, which utilize the similarity in interpretation schemes of steady state thermal and gravity fields. Its main idea was presented earlier in [13], while in the current paper we provide a complete theoretical framework of the method. Following this, we construct an algorithm for its numerical solution, which is similar to the methods employed in calculations of gravity anomalies. These results are presented in the Methods section of the current work. In the Case Study section, we test it on a three-dimensional density distribution model, which was developed earlier [17]. The problem domain extends to 80 km, which is the first regional level of the isostatic compensation. Using statistical correlation of density versus heat production and density versus thermal conductivity, we build an earth's crust and upper mantle model to solve the forward geothermal problem for density distributions given on a three-dimensional grid. Contribution of a deep thermal flux is evaluated by the method of analytical continuation to the Moho discontinuity.

Methods
The boundary value problem for the stationary heat conductivity in piecewise-homogeneous models is defined by second order differential equations. They have a strong break when the required constrains on the continuity of first derivatives are not satisfied in a space of solutions [14,16]. The contact boundaries of the neighbor media are the surfaces where the thermal conductivity coefficient has a break and, thus, a temperature gradient breaks too. Their functions have no value at these surfaces; however, the left side and right side limits exist. Such problems belong to a class of linear conjugacy with the implicitly defined boundary values (boundary conditions of the fourth kind), which are equal to the conditions of continuity of temperature and normal components of the heat flux in a 2-sided neighborhood of an ideal thermal contact [12].
In the space of generalized functions (i.e., in the linear space of continuous functionals), it is possible to avoid the construction of partial solutions for every homogeneous area and its subsequent stitching on the boundaries. A singular generalized continuous operator for an arbitrary point of the inhomogeneous space can be obtained both for the areas of continuity and the break boundaries.
The method for constructing an analytical solution to the conjugation problem in the space of generalized functions was briefly described in [13]. In Sections 2.1-2.4 (with additions in Appendices A and B), we present it in a complete form.

Generalized Continuous Operator of the Conjugacy Problem
The models of stationary thermal fields include a distribution of two parameters: the strength of radiogenic heat generation sources Q, and the thermal conductivity λ. They are connected with Geosciences 2020, 10, 199 4 of 24 an energy balance condition, namely summary flow through the closed surface S, which limits the area D τ , is equal to a summary strength of the heat sources inside this area [12]: or the same in a differential form is as follows: For the non-homogeneously heated body, in which the temperature gradient is not large, the heat flow is proportional to the first degree of the temperature gradient [12]. This is Fourier's law of heat conduction, which characterizes a conductive (or molecular) mechanism of heat transfer in a stationary media from areas with higher temperature to areas with lower temperatures. For isotropic bodies, the heat flow vector q points oppositely to the temperature gradient vector as follows: The T temperature refers to the "spreading" potential of the temperature gradients vector field (or the heat flux field). If the thermal conductivity λ is constant, the formula for the temperature matches the representation for the gravity potential with the equivalent mass density 4πQ/λ [14,15]. If the thermal conductivity is not constant, then the refracted component is added to the field of sources. At last, if the λ has an abrupt change on some boundary, then the temperature gradient component will change in the same direction. The flow q (more accurate, its normal component) is continuous on the boundary; however, every function in production (λ∇T) has a first order continuity break. In the neighborhood of such discontinuities, the boundary derivatives are replaced by the integral conditions of thermal conjugation [12].
A correctness of the continuous formulation of the conjugacy problem is ensured by a step approximation of the inverse thermal conductivity function in the area of the break surface. We can expand the operator definition area in Equations (2) and (3) by expressing its coefficients in terms of a valid combination of the inverse thermal conductivity break functions 1/λ and its gradient ∇ 1 λ as follows: Here the production of one break function (the inverse thermal conductivity and its gradient) and the continuous density function of sources Q and the normal component of the flux q exist as a valid functional combination [16,18]. Equation (2) is true not only at the points of continuity, but also at the inverse thermal conductivity coefficient break points. It corresponds to the continuous formulation of the conjugacy problem without boundary conditions. Consider the layered media D with a flat top boundary S 0 at the earth's surface and a flat bottom boundary S H (Figure 1). Inner boundaries S k (k = 1 . . . M) separate D to (M + 1) curvilinear layers D k . Above the S k boundary, the thermal conductivity of the layer is constant and equal to λ k ; the density of sources Q k is continuous or at least integrable.
Let (X) = (x, y, z) be the Cartesian coordinate system. The vertical axis (oz) is directed downwards, and the (xoy) plane acts as the top boundary S 0 . The inner boundaries S k are defined implicitly by the function Φ k (X) = z − z k (x, y) = 0. Let N k denote a unit normal vector for the surface S k . The vector N k is proportional to the gradient Φ k , and downward normal direction is selected ( Figure 1). Let ( ) = ( , , ) be the Cartesian coordinate system. The vertical axis ( ) is directed downwards, and the ( ) plane acts as the top boundary . The inner boundaries are defined implicitly by the function Φ ( ) = − ( , ) = 0. Let denote a unit normal vector for the surface . The vector is proportional to the gradient Φ , and downward normal direction is selected ( Figure 1).
Now we introduce Heaviside step functions (Φ ) . The linear combination of these step functions forms the heuristic relation of a discontinuous coefficient of the inverse thermal conductivity versus coordinates of the media.
The generalized derivative of (Φ ) is equal to the Dirac delta function (Φ ) , which is supported on the surface , and, thus, Φ ( ) = 0 [18]. Therefore, The gradient member ∇ • in Equation (4) can be converted into a form of simple layer sources as follows: A surface density of the sources ( ) is equal to the jump of normal components of the temperature gradient.
Since the heat flow is continuous, Now we introduce Heaviside step functions H(Φ k ). The linear combination of these step functions forms the heuristic relation of a discontinuous coefficient of the inverse thermal conductivity versus coordinates of the media.
The generalized derivative of H(Φ k ) is equal to the Dirac delta function δ(Φ k ), which is supported on the surface S k , and, thus, Φ k (S k ) = 0 [18]. Therefore, The gradient member ∇ 1 λ ·q in Equation (4) can be converted into a form of simple layer sources as follows: A surface density of the sources ν k (S k ) is equal to the jump of normal components of the temperature gradient.
Since the heat flow is continuous, Here, (N k ·∇T) X ∈ S k is a value of the temperature derivative, which is equal to the half sum of limit values from both sides; ε k = (λ k − λ k+1 )/(λ k + λ k+1 ) is the thermal conductivity contrast between k-th and (k + 1)-th neighbor layers. The transformed differential operator in Equation (7) is defined in break points of the thermal conductivity coefficient. Thus, as continuous functional in the generalized integral images space, it corresponds to a continual definition of the thermal coupling problem without boundary conditions.

Integral Transforms and Green's Formula
Merging of a finite number (M + 1) of layers with different thermal conductivity values forms the geothermal model of layered media with flat boundaries S 0 and S H , as shown in Figure 1.
Since the thermal coupling condition on M internal curvilinear boundaries is already included in the generalized continuous operator of the differential Equation (7), the solution should satisfy the bounding conditions only on the outer boundaries of the media. The stationary nature of the model is defined by the mixed Dirichlet-Neumann conditions [12,15]. The upper Dirichlet condition on the temperature corresponds to a isothermal mode on S 0 , while the lower Neumann condition on S H indirectly takes into account the heat flow on the boundary between the crust and the mantle.
Here θ(S 0 ) = θ 0 is the constant temperature on the upper boundary of the medium; µ(S H ) is the heat field strength at the lower boundary. It is assumed that the top and the bottom boundaries cross in infinitely remote points, making the area closed as follows: Integral image spaces of the boundary value problem (Equations (7)-(9)) is created by the convolution of the temperature Laplacian ∇ 2 T(C) and the Green's function G(A, C), which is a singular function of two variables A, C ∈ D of a homogeneous boundary value problem [12,14] as follows: Such a transform of the implicit kernel convolution (Equation (10)) is included in Green's second identity as follows: Note that the volume integral with the surface delta functions δ(Φ k ) can be transformed into a sum of two-dimensional integral forms [18] of all the internal boundaries S k as follows: Now we select the components from Green's formula, which correspond to the two types of heat flow field sources. Contribution from the internal (i.e., located inside the volume) sources Q(C) and the external sources µ(P H ) (i.e., located on the lowest plane S H ) can be combined into a sum of the temperatures, which are not distorted by the thermal conductivity contrast.
These temperatures may be considered as primary potentials. The surface simple layer sources ν k , which are created in the primary field on the inner boundaries S k , form the secondary potential: The kernel of the integral operator defined by Equations (11) and (12) is Green's function G, a known solution of a simpler boundary value problem (Equation (10)) under Dirichlet and Neumann homogeneous conditions [15]. This ensures a uniqueness of the solution if the energy balance (Equation (1)) is kept for the layered media. The difference between the integral flows on top and bottom boundaries is equal to the summary power of the inner heat sources inside the closed area D, as follows:

Low Contrast Approximation
Equation (12) for the inhomogeneous boundary problem defined by Equations (7)-(9) is a convolution of the density of the simple layer sources ν k with Green's function. Such a form of equations allows one to separate mixed boundary conditions for the layered media model. Homogeneous Dirichlet and Neumann conditions on the external boundaries S 0 and S H are included into Green's function G(A, C). Conditions of the thermal conjugation on the inner boundaries S k define the surface density of the simple layer ν k (S k ). Normal surface derivatives of Equation (12) can be substituted into a system defined by Equation (8) as follows: Here, B i (S i ) is a projection of an observation point A by N i normal to all boundaries S i ; i = 1 . . . M (see Figure 1).
Let low contrast media be defined as the media with a jump in the thermal conductivity on each S i less than 2 times (λ k+1 /λ k < 2 or λ k /λ k+1 < 2). In this case the ε i in Equation (8) is small (|ε i | < 1/3), and the integral equations system can be replaced with a linear approximation as follows: The contributions from internal and external thermal sources into the primary potential field (Equation (11) are separated. In the secondary potential they are separated as well. Substitution of the density value (Equation (15)) to the integral formula of the intermediate solution (Equation (12)) for the temperature leads us to the following formula: Geosciences 2020, 10, 199

of 24
This substitution does not change the additive structure of the operator (16). However, the kernel K(A, C) was changed as follows: The superposition of Green's functions (Equation (17)) takes into account not only the homogeneous Dirichlet and Neumann conditions, but also the heat flow refraction on its internal boundaries. Thus, the kernel K(A, C) is a full analogue for Green's function for the layered media and the flat external boundaries.

Mantle Heat Flow
The uprising component of the heat flow at the S i boundary points opposite to the z axis, which means that its sign matches the sign of the temperature gradient.
On the earth's surface S 0 (B i = B 0 ) the heat flow can be seen as a sum of the crustal radiogenic heat sources Q and the sources µ on the roof of the mantle layer q( Considering q Q (B 0 ) + q µ (B 0 ) = q observ , we obtain an integral equation for the mantle component of the heat flow.
This integral equation allows one to solve the problem of a field continuation through the inhomogeneous (by thermal conductivity) media from the earth's surface level z = 0 to the depth z = H. For the upward continuation of the mantle field component to the curvilinear boundary z = z M (x, y) we have Formulas (21) and (22) allow one to exclude the thermal field strength µ(B H ) from the foot of the mantle layer. After this, the heat flow q µ (B M ) on the roof of the upper mantle can be represented through its value q µ (B 0 ) at the earth's surface level.

Initial Data
Now we demonstrate the method using a practical example. Our study area is the Urals and the neighbor regions of Russia, which are located at 60-68 • N and 48-72 • E. This is the near-Arctic zone of an important Russian geological provinces junction: the northeastern part of the East-European platform, the Timan-Pechora plate, the northern part of the Ural Folding System and the northwestern sector of Western Siberia. Seismic survey data were used for the construction of the build frame of a 3D density model [17]. These seismic (deep seismic sounding (DSS)) velocity profiles are shown in Figure 2a, while their geographic positions compared with a map of tectonic structures [19] are presented in Figure 2b.
Geosciences 2019, 9, x FOR PEER REVIEW 9 of 23 are shown in Figure 2a, while their geographic positions compared with a map of tectonic structures [19] are presented in Figure 2b. The anomalies of the gravity field were reduced to the quasigeoid surface on a 2' × 2' grid [21,22]. A near-surface thermal flux was calculated from the temperature and thermal conductivity measurements in deep trial wells on a sparse network. This is the observed (measured) flux at the earth's surface. At this resolution the unique features of the region, such as the Urals geosyncline [23], may be seen. These features are invisible in known global databases of the observed heat flux values [24,25], which have a resolution no more than 2° × 2°. Additionally we used geothermal survey data from the Urals-Siberian region [8,26], which were interpolated from the resolution 0.5° × 0.2° (lon × lat). The data on the eastern edge of the European platform and the Timan-Pechora plate were refined using earlier publications [7,9,11,27,28]. The resulting map is provided with the current paper in the Supplementary Materials. Data of gravity and thermal measurements were recalculated to a regular grid in zone 11 of a Gauss-Kruger projection [29,30]. Figure 3a presents the anomaly gravity field for our study area on a dense 1 × 1 km 2 grid, which we also use for our density model. It is compared with the thermal flux on a sparse network of 10 × 10 km 2 , as seen in Figure 3b. The anomalies of the gravity field were reduced to the quasigeoid surface on a 2 × 2 grid [21,22]. A near-surface thermal flux was calculated from the temperature and thermal conductivity measurements in deep trial wells on a sparse network. This is the observed (measured) flux at the earth's surface. At this resolution the unique features of the region, such as the Urals geosyncline [23], may be seen. These features are invisible in known global databases of the observed heat flux values [24,25], which have a resolution no more than 2 • × 2 • . Additionally we used geothermal survey data from the Urals-Siberian region [8,26], which were interpolated from the resolution 0.5 • × 0.2 • (lon × lat). The data on the eastern edge of the European platform and the Timan-Pechora plate were refined using earlier publications [7,9,11,27,28]. The resulting map is provided with the current paper in the Supplementary Materials. Data of gravity and thermal measurements were recalculated to a regular grid in zone 11 of a Gauss-Kruger projection [29,30]. Figure 3a presents the anomaly gravity field for our study area on a dense 1 × 1 km 2 grid, which we also use for our density model. It is compared with the thermal flux on a sparse network of 10 × 10 km 2 , as seen in Figure 3b.
The main Ural fault crosses the eastern side of the Urals at approximately 60 • E. It forms the western edge of the Urals geosyncline, which is the Hercynian tectonic fold. This zone has anomalous low heat flux values [23]. It corresponds to the regional maximum of the gravity field. The similarities in morphology are traced in other tectonic zones as well. However, there is no close connection between gravity and thermal anomalies. A correlation coefficient between these maps is only −0.31, while a correlation coefficient between the relief grid and the thermal flux grid is 0.76. lat). The data on the eastern edge of the European platform and the Timan-Pechora plate were refined using earlier publications [7,9,11,27,28]. The resulting map is provided with the current paper in the Supplementary Materials. Data of gravity and thermal measurements were recalculated to a regular grid in zone 11 of a Gauss-Kruger projection [29,30]. Figure 3a presents the anomaly gravity field for our study area on a dense 1 × 1 km 2 grid, which we also use for our density model. It is compared with the thermal flux on a sparse network of 10 × 10 km 2 , as seen in Figure 3b.

Initial Density Model
Two dimensional velocity data were obtained by profile sections from 10 deep seismic sounding (DSS) profiles. Velocity cuts down to Mohorovicic discontinuity were constructed using a method of seismic tomography. Below the Mohorovicic discontinuity, there is an upper mantle, which was constructed using an isostatic compensation method and has a block structure. The blocks were selected by residual anomalies of the gravity field and then were refined by the lithostatic pressure anomaly on the 80 km depth level. This 80 km depth is the first regional level of the isostatic compensation for the Urals [17].
Velocity data were recalculated to the densities using a known piecewise-linear regression function [31]. The gravity data can be used to connect the model density with the observed gravity anomalies. The regression coefficients were clarified during the inverse linear problem solution for all the velocity profiles of the region. Thus, we obtained a possible density analogue for the velocity model for deep structures.
Density profiles were bound to a three-dimensional grid (Figure 4a). They form a carcass of a density model. The vertical resolution of the model is 100 m, so it contains 800 layers between 0 and 80 km. The horizontal step of the discretization is 1 km by 1 km. The profile data were interpolated in each layer to fill the gaps between the profiles. The average density value σ 0 was calculated for every layer. This set of density values forms the background density function σ 0 (z) of the model, which depends on depth only. Figure 4b shows σ 0 (z) along with minimal and maximal density values in the layers. As it can be seen, the average density mostly repeats the character of the minimal and maximal values. There is a sharp rise of the density at around 5-8 km below the earth's surface (the sediment cover); at middle and lower crust depths, the gradient of the density change is low, and at 40-50 km there is a clearly visible step where the mantle blocks start to appear.
layer. This set of density values forms the background density function ( ) of the model, depends on depth only. Figure 4b shows ( ) along with minimal and maximal density s in the layers. As it can be seen, the average density mostly repeats the character of the minimal aximal values. There is a sharp rise of the density at around 5-8 km below the earth's surface ediment cover); at middle and lower crust depths, the gradient of the density change is low, and 50 km there is a clearly visible step where the mantle blocks start to appear.

Refined Density Model
The method of gravity inversion was developed by the authors and was described in [32]. The inversion was performed in layers consequentially; during this process, the model was refined to satisfy the observed field. The constructed 1D background density σ 0 (z) was used for the mass continuation outside the study volume. An absolute density value in every grid element was replaced with the density jump (excess density), which is relative to the background density. The iterative scheme of corrective additions calculations in the horizontal layers provides uniqueness of the inverse problem solution [33]. The resulting density model is presented in Figure 5. It consists of 1336 × 969 × 800 prismatic elements.

Refined Density Model
The method of gravity inversion was developed by the authors and was described in [32]. The inversion was performed in layers consequentially; during this process, the model was refined to satisfy the observed field. The constructed 1D background density ( ) was used for the mass continuation outside the study volume. An absolute density value in every grid element was replaced with the density jump (excess density), which is relative to the background density. The iterative scheme of corrective additions calculations in the horizontal layers provides uniqueness of the inverse problem solution [33]. The resulting density model is presented in Figure 5. It consists of 1336 × 969 × 800 prismatic elements. The comparison of the initial distribution of minimum and maximum density values versus depth (Figure 4b) with the resulting one (Figure 5b) shows that there were almost no corrections made to the model in the lower crust and upper mantle. Therefore, the morphology of deeper boundaries (e.g., Moho), which was estimated by seismic data [3,34], cannot be refined in the process of gravity modeling. The comparison of the initial distribution of minimum and maximum density values versus depth (Figure 4b) with the resulting one (Figure 5b) shows that there were almost no corrections made to the model in the lower crust and upper mantle. Therefore, the morphology of deeper boundaries (e.g., Moho), which was estimated by seismic data [3,34], cannot be refined in the process of gravity modeling.

Grid Approximation of Thermal Physical Parameters of Geothermal Model
Petrophysics surveys evaluate the dependency between rock densities and their heat production and thermal conductivity [5,6,8,9,11]. For igneous rocks, the dependency of the density values versus heat production is most visible and was estimated by Rybach [35,36]. For sedimentary and metamorphic rocks, the dependency for the concrete region should be used [2,8,27]. Usually it is in the form of heat production variation intervals for specified density values (e.g., [5,7]). Bulashevich performed a generalization of these dependencies for the Urals region [23] (Figure 6b). It is satisfactory because for acidic to ultra-alkaline rocks, the minimum and maximum values of the heat production variation intervals are close. Using this density versus heat production dependency, the 3D density model was recalculated into a corresponding thermal model of inner sources of the radiogenic heat production field. Figure 6a presents this 3D distribution of the crustal field sources. Another important parameter of stationary thermal models is the thermal conductivity λ. Unfortunately, there is no stable and clearly visible relation between the thermal conductivity and the rock density. This happens because, in the first place, thermal conductivity depends on temperature and pressure, which partially compensate each another [23]. In addition, it depends on structural tectonic factors such as a geological age and type of stratification [5,6,9,11,37]. At this stage of research, we used values of thermal conductivity taken from a collection of rocks examples for the Urals and Siberian regions obtained by Duchkov et al. [26] and Golovanova [8]. This dependency, which takes the structural tectonic factors into account, is shown in Figure 7b. Another important parameter of stationary thermal models is the thermal conductivity λ. Unfortunately, there is no stable and clearly visible relation between the thermal conductivity and the rock density. This happens because, in the first place, thermal conductivity depends on temperature and pressure, which partially compensate each another [23]. In addition, it depends on structural tectonic factors such as a geological age and type of stratification [5,6,9,11,37]. At this stage of research, we used values of thermal conductivity taken from a collection of rocks examples for the Urals and Siberian regions obtained by Duchkov et al. [26] and Golovanova [8]. This dependency, which takes the structural tectonic factors into account, is shown in

Forward Geothermal Problem
The formulas for the stationary temperatures and heat flux for three-dimensional geothermal models of such a layered media were obtained for the condition of low contrast thermal conductivity [7,9,11]. Under this condition the calculated gravity and thermal fields are additive. This fact is the basis of our complex interpretation method for gravity and geothermal fields in inhomogeneous media.
Another important parameter of stationary thermal models is the thermal conductivity λ. Unfortunately, there is no stable and clearly visible relation between the thermal conductivity and the rock density. This happens because, in the first place, thermal conductivity depends on temperature and pressure, which partially compensate each another [23]. In addition, it depends on structural tectonic factors such as a geological age and type of stratification [5,6,9,11,37]. At this stage of research, we used values of thermal conductivity taken from a collection of rocks examples for the Urals and Siberian regions obtained by Duchkov et al. [26] and Golovanova [8]. This dependency, which takes the structural tectonic factors into account, is shown in Figure 7b.

Forward Geothermal Problem
The formulas for the stationary temperatures and heat flux for three-dimensional geothermal models of such a layered media were obtained for the condition of low contrast thermal conductivity [7,9,11]. Under this condition the calculated gravity and thermal fields are additive. This fact is the basis of our complex interpretation method for gravity and geothermal fields in inhomogeneous

Results
Three dimensional distributions of the heat production sources Q in the earth's crust and thermal conductivity λ in crust and mantle in the form of a 3D array are shown in Figures 6 and 7 correspondingly. The heat flux field was calculated from these data using Equation (19). The Q/λ component of Equation (19) was precalculated by an elementwise division of Q(x k , y k , z k ) array (Figure 6a) by λ(x k , y k , z k ) array (Figure 7a). The corrections for the thermal conductivity contrast from the simple layer sources ν(S k ) were calculated for the model of two boundaries, and, thus, three layers [34]: a sediments layer with low thermal conductivity (λ ≈ 2 w·m/K), and crustal (λ ≈ 2.6 w·m/K) and mantle layers with high thermal conductivity (λ ≈ 4.2 w·m/K). Green's function G for the two flat boundaries (Equation (A10), Appendix A) and the function K for the layered media (Equation (A17), Appendix B) correspond to the function of inverse distances in the gravity problem. Figure 8 presents the solution of the forward thermal flux calculation for the model using Equations (19) and (20). The heat flux for the earth's crust is calculated at the earth's surface level z = 0. The mantle component of the flow (Equation (20)) is calculated from the difference between the observed and calculated values of the heat flow.
The mantle component of the flux at the upper boundary of the mantle is changing in wide interval from the normal (average) values of 20-30 mW/m 2 at shields and ancient platforms (e.g., Eastern European platform and Timan Pechora plate) [2,5,27]; to high values of 40-50 mW/m 2 and more at the western end of the young Western Siberian plate [11,37]. The mantle flux at the Urals is partially negative. On the one hand this is connected to the extremely low values of the measured temperature gradients and, thus, the calculated heat flux, in wells of the Urals geosyncline [8,23]. On the other hand, the relatively high thermal production at the mountain roots of the Urals Folding System may be the reason [9]. The negative values of stationary thermal flux at the earth's surface level are against the physics of geothermal processes and require a reasonable adjustment of measured data taking the paleoclimate into account [38,39]. function for the two flat boundaries (Equation (A10), Appendix A) and the function for the layered media (Equation (B7), Appendix B) correspond to the function of inverse distances in the gravity problem. Figure 8 presents the solution of the forward thermal flux calculation for the model using Equations (19) and (20). The heat flux for the earth's crust is calculated at the earth's surface level = 0. The mantle component of the flow (Equation (20)) is calculated from the difference between the observed and calculated values of the heat flow. The mantle component of the flux at the upper boundary of the mantle is changing in wide interval from the normal (average) values of 20-30 mW/m 2 at shields and ancient platforms (e.g., Eastern European platform and Timan Pechora plate) [2,5,27]; to high values of 40-50 mW/m 2 and more at the western end of the young Western Siberian plate [11,37]. The mantle flux at the Urals is partially negative. On the one hand this is connected to the extremely low values of the measured temperature gradients and, thus, the calculated heat flux, in wells of the Urals geosyncline [8,23]. On the other hand, the relatively high thermal production at the mountain roots of the Urals Folding System may be the reason [9]. The negative values of stationary thermal flux at the earth's surface level are against the physics of geothermal processes and require a reasonable adjustment of measured data taking the paleoclimate into account [38,39].

Inverse Problem of Analytical Fields Continuation
A deep heat flux, which is given through the Neumann boundary problem (Equation (9)), is a component of the solution of the forward problem for the temperatures (Equation (16)) and heat flux (Equation (20)). However, it is not an a priori known function. Nevertheless, a "trace" of this function can be clearly seen in the field flux difference (observed minus radiogenic) at the upper boundary z = 0 (Figure 8c). The lower bounding surface z = Н is a helpful carrier of simple field sources with a density of ( , ). Thus, the right hand part of Equation (21) is a convolutional integral operator for a direct recalculation of the field from the level = to the level = 0. An inversion of the operator of Equation (21) gives us the solution of the inverse problem of the analytical continuation

Inverse Problem of Analytical Fields Continuation
A deep heat flux, which is given through the Neumann boundary problem (Equation (9)), is a component of the solution of the forward problem for the temperatures (Equation (16)) and heat flux (Equation (20)). However, it is not an a priori known function. Nevertheless, a "trace" of this function can be clearly seen in the field flux difference (observed minus radiogenic) at the upper boundary z = 0 (Figure 8c). The lower bounding surface z = H is a helpful carrier of simple field sources with a density of µ(x H , y H ). Thus, the right hand part of Equation (21) is a convolutional integral operator for a direct recalculation of the field from the level z = H to the level z = 0. An inversion of the operator of Equation (21) gives us the solution of the inverse problem of the analytical continuation of the mantle component of the heat flux down through layered inhomogeneous media in the direction of field sources.
The inverse problem of the fields continuation is unstable with respect to irregular noise in the input data. This is why the algorithm for the analytical continuation should include a regularization parameter. Its value is selected that way so the continuation of the regularized solution to the initial height level matches (with some error) the initial field. The uprising heat flux at the upper isothermal surface has a vertical component only. The heat component of earth crust sources is excluded, and the remaining mantle component is recalculated to a vertical gradient of temperature µ(x, y) at the lower bounding surface z = H. This is enough for the ability to correct geothermal parameters of boundary function and calculation of the temperature and heat flux vector components at every point of the layered inhomogeneous media.
To match the deep mantle component of heat flux with the relief of the upper mantle roof, the components of the temperature gradient vector were calculated as well as the uprising heat flux (Equation (22)) at the depths, which correspond to the Mohorovicic boundary position M (x,y). Figure 9 presents the results of the vertical component of mantle flux recalculation to the curvilinear surface of M boundary. The morphology of the recalculated fields repeats, in general, the relief of mantle structures, and they have a good match each with the other. The correlation coefficient between the mantle relief grid and thermal flux grid is 0.76, which was calculated by 133 × 97 grid elements on q observ (x, y, 0) field grid.
The increased heat flux corresponds to the uplift (an increase in elevation) of Moho on the west side of the Western-Siberian plate [37], while the low flux values correspond to the boundary depression along the Urals mountain belt (approximately 59 • E) and neighbor foreland basins [9]. On the contrary to these anomalous geothermal regions, the mantle flux at the eastern edge of the Eastern European platform and northern part of the Timan-Pechora plate is close to its background average values. This is usual for ancient shields and platforms [7,11,27].
The heat flux and the vertical component of the calculated temperature gradient at Moho (Figure 9b) correlates with the age of the geological structures of the earth's crust as well as with the seismic waves velocity change at the Mohorovicic discontinuity [5,7]. Such a velocity inhomogeneity must be accompanied with inhomogeneous temperature distribution [10,36]. The lower the seismic waves velocity along the mantle boundary M is, the higher the temperature values of the mantle blocks [2,3,11,37].
boundary function and calculation of the temperature and heat flux vector components at every point of the layered inhomogeneous media.
To match the deep mantle component of heat flux with the relief of the upper mantle roof, the components of the temperature gradient vector were calculated as well as the uprising heat flux (Equation (22)) at the depths, which correspond to the Mohorovicic boundary position M (x,y). Figure  9 presents the results of the vertical component of mantle flux recalculation to the curvilinear surface of M boundary. The morphology of the recalculated fields repeats, in general, the relief of mantle structures, and they have a good match each with the other. The correlation coefficient between the mantle relief grid and thermal flux grid is 0.76, which was calculated by 133 × 97 grid elements on ( , , 0) field grid. The increased heat flux corresponds to the uplift (an increase in elevation) of Moho on the west side of the Western-Siberian plate [37], while the low flux values correspond to the boundary depression along the Urals mountain belt (approximately 59° E) and neighbor foreland basins [9]. On the contrary to these anomalous geothermal regions, the mantle flux at the eastern edge of the Eastern European platform and northern part of the Timan-Pechora plate is close to its background average values. This is usual for ancient shields and platforms [7,11,27].
The heat flux and the vertical component of the calculated temperature gradient at Moho ( Figure  9b) correlates with the age of the geological structures of the earth's crust as well as with the seismic waves velocity change at the Mohorovicic discontinuity [5,7]. Such a velocity inhomogeneity must be accompanied with inhomogeneous temperature distribution [10,36]. The lower the seismic waves velocity along the mantle boundary M is, the higher the temperature values of the mantle blocks [2,3,11,37].
An isostasy, which has been taken into account in the mathematic model, the change of the density value is controlled. The blocks in the model of the mantle in Figure 4 were selected taking the isostatic compensation hypothesis into account. The density of mantle blocks depends on the Moho boundary relief. Uplifts in the Moho correspond to the increase in thickness of the isostatic An isostasy, which has been taken into account in the mathematic model, the change of the density value is controlled. The blocks in the model of the mantle in Figure 4 were selected taking the isostatic compensation hypothesis into account. The density of mantle blocks depends on the Moho boundary relief. Uplifts in the Moho correspond to the increase in thickness of the isostatic compensation layer and the decrease of the density (and the velocity) in a mantle block. The converse is also true; a decreased Moho relief corresponds to increased values of density and velocity in underlying mantle blocks [17].
All the required elements for the evaluation of the mantle heat flux were clarified by authors earlier during the construction of the velocity and density models for the region [17,[31][32][33][34]. The inner agreement of the results obtained by gravity and thermal fields does not contradict with the deep seismic sounding data. Thus, this presented method is a simple and quite informative technique for the joint interpretation of geothermal, gravity and seismic data for three dimensional mathematical modeling for the earth's lower crust and upper mantle studies.

Conclusions
The theory of thermal conjugacy problems solution is based on an analytical approximation of the inverse thermal conductivity coefficient, which has breaks, at the thermal contact surfaces. A differential operator for the equation of stationary thermal conductivity in the space of generalized functions is continuous. A transformation of this continuous operator convolution has given us the solution of the thermal conjugacy problem in the form of a single integral formula.
An analytical form solution for the Dirichlet-Neumann mixed edge value problem was obtained for the model of the layer with thermal conductivity jumps on boundaries in the case of an integrable function of radiogenic heat flow sources distribution. At the break surfaces of the thermal conductivity, the equivalent sources of the simple layer are formed. Their densities are found from the Fredholm integral equation system of second order. For the media, which have low contrast in thermal conductivity, a density of simple layer sources is found in the analytical form and the solution for the heat flow becomes additive.
We have developed and implemented the algorithms for forward and inverse problems of stationary thermal conduction for layered media. The principle of field superposition is the basis for interpretation of thermal and gravity anomalies.
Using these algorithms, the modeling of thermal parameter distribution in the earth's crust was performed. The forward thermal problem solution was obtained. The heat flux and temperature calculations were made for the model of the density distribution of the Ural region's upper lithosphere.
The mantle component of the heat flux at the earth's surface was estimated. An algorithm for the continuation of this mantle component of the heat flux from the earth's surface level down to the Mohorovicic discontinuity surface was developed. A tectonic zoning was performed for the obtained map of the mantle component of heat flow. The relief of the Moho boundary was compared to the mantle heat flux anomalies.
Thus, the correctly chosen strategy of a continuum approach to the boundary problem conjugation for the models of piecewise-homogeneous media (the edge value problems with conditions of the fourth order) allows one to obtain the solution for quite a complex media geometry in the form of explicit mathematical expression. At the end, it allowed us to construct the unified additive algorithm for the stationary thermal field calculation in inhomogeneous media. Joint usage of gravity and thermal field data in the interpretation improves the reliability of the modeling and reduces the degree of non-uniqueness of the results of gravity and geothermal data inversion.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Green's function G is a solution for the Poisson equation with singularity in the right hand part as follows: Here, X = (x, y, z) is the Cartesian coordinates of current parameter point, and X 0 = (x 0 , y 0 , z 0 ) is the Cartesian coordinates of source point. The solution is sought in the horizontal strip X ∈ (−∞; +∞) × (−∞; +∞) × [0; H] for the inner localization of the source X 0 = (x 0 , y 0 ; 0 < z 0 < H). Homogeneous boundary conditions corresponds to the mixed Dirichlet-Neumann boundary value problem as follows: Let X = (ρ, φ, z); X 0 = (ρ 0 , φ 0 , z 0 ) be the auxiliary cylindrical coordinate system with depth axis z located at the same position as depth axis of the source.
The density of the singular source in Equation (A1) can be transformed into an axially symmetric function as follows: With these new variables (ρ − ρ 0 ) = η; z = z, the problem has an axial symmetry, and the G function does not depend on the angle φ. Therefore, Let F(s, z) be the Hankel transform of G(η, z) with Bessel functions for order 0 [40] as follows: In the space of the integral transformants, Equation (A3) is replaced with the ordinary differential equation of the second order (J 0 (0) = 1) : In the space of images, the boundary conditions (Equation (A2)) are changed as The solution of inhomogeneous Equation (A5) is a sum of a complementary solution and the partial solution, V(s, z), of Equation (A4) as follows: where A(s) and B(s) are free coefficients, which depend on s, and V(s, z) is a Duhamel integral Using boundary conditions (Equation (A5)), we can obtain A(s) and B(s) coefficients The duality of the partial solution (Equation (A6)) leads to the duality of the analytical representation of a general solution for images as follows: .
Using obvious transformations, hyperbolic functions are converted to exponential series as follows: where z < = min(z, z 0 ); z > = max(z, z 0 ); A preimage of the function, G(η, z), is calculated using the inverse Hankel transformation [40] as follows: Each member of the exponential series of images (Equation (A8)) is converted into a series of corresponding preimages for the inverse distance functions [41] as follows: Alternating series for Green's function for the mixed Dirichlet-Neumann problem can be interpreted as a series of sequential mirror (or electrostatic) reflections [14]. For the odd reflections of a source from the isothermal (or conductive) plane, the sign is alternating. The even reflections are from the thermally isolated (or non-conductive) plane and in this case the sign is kept. Thus, for every source inside the layer z 0 ∈ [0; H], the (4n + 1) images of the sources are formed, which are located in upper (z < 0) and lower (z > H > 0) semispaces. These images are represented with series in Equation (A9). If we extend the n index to the negative values, the solution formula (Equation (A9)) will have a universal and compact form as follows: Here, R is a distance between a source point and the observation point in projection to the horizontal plane: R = ρ − ρ 0 = (x − x 0 ) 2 + (y − y 0 ) 2 . It seems that in the form of Equation (A10), Green's function was initially obtained in the electric survey problem when a sedimentary structure of colorful Nobili rings was studied [15]. The boundary conditions of these electric and geothermal problems are equal.

Appendix B
Consider a thermal conjugacy problem for the layered media with inhomogeneous Dirichlet-Neumann conditions at the external boundaries. Its solution is an integral formula (Equation (16)), which is a convolution of the internal (i.e., distributed in the volume), Q(τ C ), and the external (i.e., located on the surface), µ(S H ), field sources with Green's function K(A, C). K(A, C) is a solution (Equation (17)) of the same conjugacy problem, but for the singular density of sources.
Equation (A11) in this form separates mixed boundary conditions for the layered model. The Dirichlet-Neumann conditions at the external boundaries S 0 and S H are used in Green's function G(A, C), while the conditions of the thermal conjugacy at the thermal contact surfaces S k are implemented by convolution of Green's function and its normal derivative at the internal boundaries.
At the isothermal boundary z = 0 (earth's surface), an uprising heat flux (Equations (19) and (20)) is proportional to the vertical component of the temperature gradient. The latter is calculated as a vertical gradient of Green's function K at point A.
Here, ε k is the thermal conductivity contrast at S k .
Using the Cartesian coordinate system (X) = (x, y, z), we define the distances between a source (or its image) point, C, and two points, namely an observation point, A, and a contact surface S k point, P k = P.
Here, R AC , R PC are the distances between corresponding points in the horizontal projection: Let us rewrite a series (Equation (A10)) of Green's function G(A, C) using this new notation.
A vertical derivative G(A, C) with respect to A point (n → m) yields the following: A similar method is used for the calculation of the Green's function G(P k , C) gradient with respect to the observation point P k = P. However, at the curvilinear surface S k , there are all three components of the gradient vector, so the convolution of normal Green's function derivatives can be, in the general case, found numerically only.
Thermal conductivity altering between layers causes effects of the thermal shielding. Its value can be estimated by a simplified model of horizontal layers. Let S k be the horizontal plane z = z k .
Geosciences 2020, 10, 199 20 of 24 A normal to this plane is directed towards the depth axis. The directional derivative of G with respect to this normal has the vertical component only: (N P ·∇ P )G(P, C) = ∂G(P, C) ∂z P = − ∂ ∂z C +∞ n = −∞ (−1) n 1 (X PC (z P + 2nH − z C )) + 1 X PC (z P + 2nH + z C ) (A14) For the flat boundaries S k the integral convolution (Equation (A12)) is calculated analytically. Now we introduce the auxiliary cylindrical coordinate system, whose vertical axis goes through the point A.
To calculate the inner integral of Equation (A16) with respect to the angle φ, we use a formula for the Bessel functions with integer order [41] as follows: Here, sgn( η m ) is a symmetric function of its variable, η m , sign: The result of the convolution calculation in Equation (A12) for the derivative K z can be compactly written using symbolic from Equation (A15) as follows: The substitution of η m = η ∓ m values results in the sum of [(2m + 1) + (2m + 1)2(2n + 1)] = (2m + 1)(4n + 3) series elements. δ m here means that the element at m = 0 is taken with the coefficient 1/2: The series (Equation (A17)) for the K z (m, n) derivative is a kernel for the integral operator of forward problem (Equations (19) and (20)) and the problem of analytical field continuation downwards to the depth H (Equation (21)). If z P ≤ H/2 then we can ignore the lower boundary effect and get a simpler formula by setting m = n = 0 in Equation (A17) as follows: 2z Pk +z C [R 2 AC +(2z Pk +z C ) 2 ] 3/2 − sgn(z Pk − z C ) (z Pk +|z Pk −z C|) R 2 AC +(z Pk +|z Pk −z C|) From Equation (A18), it is clearly seen that z C sources below and above the contact surface z Pk behave differently in the field of temperature gradients K z . For example, an income from µ(z H ) component into the subsurface heat flow (Equation (20)) is calculated from the deep level only, i.e., z H = z C > z Pk . In this case an effect of the thermal shielding is explicit: Depending on the ε k sign, the temperature gradient on the isothermal plane is higher (for negative) or lower (for positive) than that of the homogeneous media. The same scheme (Equation (A18)) describes the shielding of the thermal field from sources z C by upper thermal contact boundaries z Pk .