Multiscale Multiphysics Modeling of the Inﬁltration Process in the Permafrost

: In this work, we design a multiscale simulation method based on the Generalized Multiscale Finite Element Method (GMsFEM) for numerical modeling of ﬂuid seepage under permafrost condition in heterogeneous soils. The complex multiphysical model consists of the coupled Richards equation and the Stefan problem. These problems often contain heterogeneities due to variations of soil properties. For this reason, we design coarse-grid spaces for the multiphysical problem and design special algorithms for solving the overall problem. A numerical method has been tested on two- and three-dimensional model problems. A a quasi-real geometry with a complex surface is considered for the three-dimensional case. We demonstrate the efﬁciency and accuracy of the proposed method using several representative numerical results.


Introduction
It is very important to consider the permafrost strata in soils under difficult climatic conditions in the Far North and the Arctic. Moreover, the structure is not constant, as many may believe, but is subject to changes as a result of thawing and freezing processes. Critical changes in structure, as a result, can lead to catastrophic consequences. Furthermore, the case with the fuel storage company NorNickel clearly demonstrates it. According to the official version, an accident occurred due to a strong thawing of the permafrost layer that led to the damage as the base could not withstand the weight of the structure. It caused irreparable ecological damage in the region and incurred colossal financial losses. In the permafrost zone, one of the main tasks in construction of buildings is a reasonable determination of the bearing capacity of permafrost soils. In addition, from a technological point of view, it is crucial to plan the start period of construction. Construction start date depends on the results of geocryological forecast and thermal engineering calculation taking into account many factors such as structure, composition, and soil properties, meteorological conditions, temperature of each engineering and geological element, the thermal effect of buildings and structures during operation. When carrying out the predictive calculations for construction under permafrost conditions, one of the important tasks is to determine the three-dimensional thermal state of soils during the operation of the structure, since the stability and reliability of the structure depend on the depth of thawing and soil temperatures. The climate is also one of the main factors in the dynamics of physical processes.
The model equations can be generalized and include sink terms due to, for example, root uptake [1][2][3]. In this case, one needs new multiscale basis functions that describe root uptake and this has not been studied in the literature. Our goal in this paper is to focus on construction related issues; however, the additional effects due to water uptake in some applications. In the paper, we design and implement Generalized Multiscale Finite Element Method (GMsFEM) for flows into heterogeneous permafrost soils. We construct a mathematical model by combining several models [4][5][6]. The seepage process is implemented using the Richards Equation [7][8][9][10][11][12][13][14][15], where the coefficients of permeability and the derivative of saturation concerning pressure are empirical dependences based on estimates of the percolation rate. The process of heat transfer in the soil is described by the heat conduction equation and takes into account the phase transition of pore moisture into ice and vice versa [16,17]. An additional convective term introduced considers the effect of saturation on temperature, and the effect of temperature on the seepage process is taken into account through the permeability coefficient [5].
Due to the multiscale nature of the problem, direct numerical simulations can be resource-intensive. For this reason, we introduce some types of upscaling or multiscale methods. These methods solve the problem on a coarse grid by introducing effective media properties or multiscale basis functions (e.g., [18,19]). The extensions of these methods to complex multiphysical problems require some special treatments. In this paper, we design GMsFEM methods for our coupled multiphysical problems. The computational algorithm is based on the GMsFEM [18][19][20][21][22][23].
We would like to highlight our contributions. In [24,25], we have considered a linear heat transfer, where the permeability does not depend on the pressure. This simplifies the algorithm as one does not need to iterate and update multiscale basis functions. In this paper, we consider the nonlinear soil model, which is more realistic. In this case, the permeability depends on the pressure and the overall multiscale procedure requires a somewhat different approach. We design multiscale basis functions and iterative methods for solving the global multiphysical problem.
Multiscale methods have become very popular in recent years and a variety of methods was developed, for example, Multiscale Finite Volume Method, Heterogeneous Multiscale Method, Multiscale Finite Element Method, Variational Multiscale Method, and so on [26][27][28][29]. For high contrast porous media, more than one degree of freedom should be introduced for accurate approximation of the processes. The main idea of GMsFEM is to apply multiscale basis functions to obtain important information in each coarse grid (computational grid) and build a reduced-order model. In this method, we construct a coarse mesh then compute the snapshot space at each coarse element and construct multiscale basis functions by performing the appropriate local spectral decomposition in each coarse block. The types of local spectral problems are motivated by analysis.
In the paper, we present several two-dimensional and three-dimensional numerical tests. In our tests, we choose parameters and test simulations by using different number of basis functions per each coarse-grid block. Our results show that using fewer basis functions, one can achieve a reasonably accurate approximation of the solution.
The work consists of 5 chapters and an introduction. The second chapter contains the statement of the problem. It discusses the process of water seepage into frozen ground. The third chapter provides a finite element approximation of the calculated mathematical model. In the fourth chapter, we demonstrate GMsFEM. The last two chapters provide numerical results for a 2D and 3D problem. The paper ends with the conclusions based on the results of calculations.

Mathematical Model
We consider the process of water infiltration into the ground under permafrost conditions. To do this we write down the associated mathematical model: Seepage process. To describe the seepage process we use the Richards equation that generalizes Darcy's law. Note that there are three different forms of writing the Richards Equation [9,10]: in terms of pressure, in terms of saturation, and mixed form. We in turn use the Richards equation written in terms of pressure: here, p =p/ρg is head pressure,p is pressure, m is porosity, s(p) is saturation, K(p) is hydraulic conductivity.
The following dependencies are true for the coefficients: where K s is fully saturated conductivity, γ, σ are problem coefficients.
Heat transfer process. To simulate the thermal regime of soils, we consider which thermal conductivity equation is used, taking into account the phase transitions of pore moisture. In practice, phase transformations occur in a small temperature range [T * − ∆, T * + ∆]. Let us take sufficiently smooth functions φ ∆ and δ ∆ (T − T * ) depending on temperature: Then, we obtain the following equation for the temperature in the region Ω: and L is specific heat of phase transition (the latent heat). The resulting Equation (4) is a standard quasilinear parabolic equation.
For the coefficients of the equation, the following relations are true here, ρ + , c + , λ + , ρ − , c − , λ − are density, specific heat, thermal conductivity of thawed and frozen zones, respectively. Fully coupled. We adapt the complete physical model by analogy with [5]. The effect of saturation on temperature is taken into account by introducing an additional convective term: The effect of temperature on the seepage process is taken into account through the permeability coefficient (if we mark the hydraulic permeability through K(p)): here, ε = 10 −6 is small number. Thus, based on (1), (2), (4), (6), (7), we write down the complete system of equations describing the seepage process in a porous medium, taking into account temperature and phase transitions.
Boundary and initial conditions. We consider a quasi-real domain Ω ⊂ R 2 , with boundary Γ = ∂Ω, Figure 1). Let us supplement the complete system with boundary and initial conditions: For temperature. On top of the area (Γ st ∪ Γ in ): here, T air is air temperature.
On the lateral and lower parts of the border: For pressure (saturation) there are non-flow conditions everywhere: except part Γ in , where under the condition temperature of the air is greater than fifteen degrees of Celsius: Initial conditions for temperature: Initial conditions for pressure (saturation):

Fine Grid Finite Element Approximation and Picard Iteration for Linearization
To describe the approximation in time we introduce a uniform grid with a step τ: and introduce the notation p n = p(t n ), T n = T(t n ). For approximation in time, we will use an analog of the implicit difference scheme. For the Richards equation, we use the Picard iterations, and for the Stefan equation, linearization from the previous time layer is used.
where k are Picard iterations. Only after solving the equation for p can we calculate the equation for T using p n+1 = p k+1 .
For the approximation in space, standard Lagrangian finite elements are used. Next, we introduce the space of finite elements: Here H 1 (Ω)-Sobolev space. Next, we represent the system in a variational form (16). Thus, we obtain a variational formulation with taking into account the boundary conditions: here, v 1 ∈ V h p and v 2 ∈ V h T .

Generalized Multiscale Finite Element Method
In this section, we describe how to build a local reduced-order model on a snapshot space by solving some local spectral problems using GMsFEM. First, we need to generate a rough mesh and build a snapshot space. Then, we solve special spectral problems in these snapshot spaces for each coarse grid block we get some kind of set of multiscale basis functions. Thus, it can be seen that the key component of the GMsFEM method is the construction of local basis functions.
In offline computation, we have to build a snapshot of the V ω i snap space for each rough neighborhood ω i . The snapshot space can be the space of all small-scale basis functions or solutions of some local problems on ω i with all possible boundary conditions.
here, δ j (x) some set of functions defined on ∂ω i , j = 1, . . . , J i (J i number of all nodes ∂ω i ). Therefore, we define Offline space is constructed using the following local spectral problems in snapshot space: Here A and S denote similar fine-scale matrices defined as follows: here, φ m are linear basis functions.
To generate the offline space, choose the smallest M ω i off eigenvalues and find the corresponding eigenvectors ψ off k = ∑ m Ψ off k,m φ off m for k = 1, . . . , M ω i off .
We implement special multiscale basis functions φ i j (Figure 2) to describe near surface form effects. We must multiply the found eigenvectors by the partition of unity functions χ i .
off denotes the number of eigenvectors that are sampled for each local area ω i . To define a partition of unity function, we first define an initial coarse space ; here, N the number of rough neighborhoods and χ i is a standard multiscale partition of unity function which is defined by: where g i is a continuous function on ∂C and linear on each edge ∂C; here, C is the cell of the coarse grid.
Next, we define a multiscale space: and define the projection matrix: In this problem, obtained basis functions are used to solve a fully coupled problem. Using the projection matrix R, we solve the problem using a coarse grid: where M c = RNR T , A c = RAR T , F c = RF and u n ms = R T u n c ; here, u n ms is a fine-grid projection of the coarse-grid solution. M and A are the mass and stiffness matrices for the Fine system, respectively, F is the vector of the right-hand side and u is the required function for the pressure P and T.

Numerical Results Two-Dimensional Problem
Numerical simulation of an applied problem in a two-dimensional formulation describing water seepage into the permafrost. The object dimension is 10 m wide and 5 m deep ( Figure 3). In an open area boundary conditions of the third kind were used-the external environment. For the parameters of the external environment, the monthly average values of air temperature were taken in the area of Yakutsk for the last year (Table 1). In these calculations, it was assumed that the heat flow from the bowels would not affect the temperature distribution of the rocks; therefore, the homogeneous Neumann boundary condition was taken at the lower boundary of the computational domain.
We implement numerical modeling of the problem under consideration for the following values of the thermophysical properties of the soil: For numerical comparison of the fine-scale and multiscale solutions, we use weighted relative L 2 and energy errors for temperature and pressure: where u h and u ms are the fine-scale and multiscale solutions. Table 2 shows the relative errors of L2 and energies for a different number of multiscale basis functions. First of all, we noticed that by updating the basis functions more often we can get more accurate solutions. We obtain 0.82% L2 error for the pressure and 1.3% L2 error for temperature on 2 multiscale basis functions. On the other hand accuracy of the method increase on 8 multiscale basis functions. In this case, method provide 0.16% L2 error for the pressure, and 0.23% for temperature. The coarse grid solution using 8 basis functions for each temperature and pressure are shown in Figures 4 and 5 for four time steps.

Numerical Results Three-Dimensional Problem
We expand the 2D problem to the problem in a three-dimensional setting ( Figure 6). The area in the plan has the same dimensions of 10 m and a height of 6 m. The computational grid has the dimensions N n = 522,774 and N e = 35,844,142. All characteristics of the problem remain the same as in the case of 2D. The calculations were carried out for 1 year with a time step of τ = 24 h. To generate the soil surface, we used the following surface equation z(x, y) = 5.5 + 0.5 sin(x + y) − 0.2 exp(−0.5[(5 − x) 2 − (5 − y) 2 ]/10). At the center of this geometry, there is a pronounced depression through which liquid seeps. This depression serves as an analog of places where water from precipitation accumulates.  Table 3 demonstrates numerical convergences in norm ||e|| L 2 and ||e|| a for temperature and pressure. Methodical results for the three-dimensional case qualitatively repeat calculations in two-dimensional calculations. At the same time, the main trends persist and a decrease in the error can be observed with an increase in the number of multiscale basis functions. The main computational difficulty is localized in the Richards equation. This fact can be described by the error in ||e|| a norm. This is due to the non-linearity of the equation complicated by the complex permeability coefficient which depends on the temperature. This can be dealt with by an increasing number of multiscale basis embassies. On the other hand, the average ||e|| L 2 error demonstrates good accuracy of the applied method. Already good accuracy is achieved when 16 multiscale basis functions are implemented. We obtain 0.52% L2 error for the pressure and 0.18% L2 error for temperature on 16 multiscale basis functions. Discussed L2 norms are observed to provide smaller errors because they do not contain gradients. It is known that the gradients of multiscale functions are rough, and these spatial fields are more difficult to represent with multiscale methods. The coarse grid solution using 8 basis functions for each temperature and pressure are shown in Figures 7 and 8 for four time steps. Multiscale solvers can significantly reduce the size of the system and provide accurate solutions.  These results indicate that our method is robust with respect to the contrast in the coefficient, and is able to give accurate approximate solution with a few local basis functions per each coarse neighborhood. Numerical results demonstrate fact that the infiltration process strongly affects the frozen ground.

Conclusions
A generalized multiscale method for solving the problem of the seepage process into permafrost soil is presented. Such kinds of problems are relevant for applied problems which involve the processes of thawing and freezing of a permafrost layer. The adaptive basis functions have been developed that take into account irregularities on discretization level for the complex geometry of the surface. The Multiphysics model was assembled based on two nonlinear problems (Richards equation and Stefan problem) for numerical implementation. We want to that our results are numerical and further studies are needed to obtain the convergence. Based on the foregoing, the proposed method has shown its efficiency both in simplified two-dimensional problems and in applied real threedimensional cases. To demonstrate modeling potential the results of numerical calculations carried out in conditions close to the Yakutia region are presented.