Filtration of Radioactive Solutions in Jointy Layers

1. Introduction

The solution of global energy problems of mankind, first and foremost associated with the development of nuclear energy. Already by 2030 the share of nuclear power generation in total electricity production should be about 25-30% (today - 16%). Currently, the total amount of radioactive waste in Russia is estimated at 5 10⁸ m³, the total β-activity of which is estimated at 7.3 10¹⁹ Bq. At the same time on the liquid radioactive waste (LRW) accounts for about 85% of total activity, and their treatment and disposal become the most important task of nuclear energy.

One of the safest ways of disposal of waste of nuclear and chemical production is injection of them into deep-seated subterranean formations. Therefore, an important issue is to study the processes of the joint heat and mass transfer during the injection of waste into a porous collector layer to predict and control the state of the areas covered by the influence of radioactive impurities. The above forecast is carried out mainly by calculations, since the possibility of experimentally sizing of deep zones of contamination is very limited.

The processes of mass transfer in porous media have long been the object of study for many researchers. Have become classics of the G.I. Barenblatt [1], Bear J. [2, 3, 4], Bachmat Y. [5], A. A. Ilyushin [6], V. M. Keyes [7], L. D. Landau [8], R. I. Nigmatullin [9], V. N. Nikolayevsky [10, 11], L. I. Sedov [12]. In the works of Prakash A. [13], A. A. Barmin [14], E. A. Bondarev [15], M. L. Zhemzhurov [16, 17], E. V. Venetsianov, R. N. Rubinstein [18] the problems of the filtration of solutions, taking into account the phenomenon of adsorption, are regarded. Fluid flow through porous materials [19-24] are coincided to be well studied. Study of models of multi-component flows is devoted to the work of R. E. Swing [25].

Problems of disposal of radioactive waste in geological formations and the resulting ecological problems discussed in works of A. S. Belitsky, E. Orlova [26], A. Rybalchenko, M. K. Pimenov [27]. Modeling of temperature and radiation fields examined in works of D. M. Noskov, A. D. Istomin, A. G. Kessler, A. Zhiganov [28-29] (Seversk Technological Institute), I. Kosareva, and E. V. Zakharova [30-31] (Institute of Physical Chemistry RAS), and other researchers. In the works of A. Lehova, Y. Shvarova studied the rate of radionuclides in groundwater, the behavior of radioactive waste in the earth's crust after the injection. At the same time remain relevant problem of determining the concentration dependence of the fields on the parameters of injection of radioactive impurities, injection technology on the parameters of layers, etc.

The study of filtration processes in multilayer formations, as well as any thermodynamic problems of contacts of the bodies and environments, leads to the necessity of solving the problems of conjugation. To solve these problems are widely used numerical methods. The analytical solutions are constructed only for simple cases, such as linear flow in mass-isolated formation [32, 33]. And as the disposal at the request of the IAEA carried out on the timing of the order of tens thousand years, in these circumstances, the porous layer, can hardly be considered mass-isolated.

In this paper, in example of study of the filtration process of radioactive solutions, represented a modification of the asymptotic method, allowing successfully construct approximated solutions to conjugacy problems.

2. The mathematical formulation of the problem of heat and mass transfer in fluid flow with radioactive contaminant in the deep layers

Let us consider problem of heat and mass transfer, which describe the interrelated fields of concentration and temperature of the radioactive contaminant in the porous layer, through which flows a liquid with impurities, and the covering and the underlying layers are waterproof.

Typically, in deep horizons an aqueous solution is injected. This solution consist of a different soluble chemical compounds formed during the acid treatment process of structural elements of reactors and other parts of the design (process waste), or in the decontamination of buildings, cars, clothing and so on (non-technological waste) and includes a mixture of various radioactive nuclides [34]. Quite naturally the initial density of the solution divided into two factions

$${\rho }_{tot}={\displaystyle \sum _{k=0}^{Nch}{\rho }_k}+{\displaystyle \sum _{i=1}^{Nrc}{\rho }_i},$$E1

where $${\rho }_k$$ is the density of the dissolved non-radioactive components (for k = 0, we obtain the density of the solvent) $${\rho }_i$$ is the concentration of radioactive i-th nuclide, N_ch, N_rc - the number of different non-radioactive and radioactive components in the solution, respectively.

Consider an arbitrary reaction volume $$dV$$ in a porous layer containing a multi-component mixture (1). Mass flow passing through the surface of dS reaction volume $$dV$$ can be represented as the sum of four terms

$$\frac{\partial }{\partial \tau }{\displaystyle \int {\rho }_{j}dV}={\displaystyle \underset{S}{\int }{\overrightarrow{j}}_j\overrightarrow{n}dS}+{\Delta }_j^{ch}+{\Delta }_j^{reac}+{\Delta }_j^{ex},$$E2

where the first term takes into account the mass exchange with the environment through the diffusion and convection currents, the second term describes the rate of change of mass in chemical reactions, the third term takes into account the change in mass due to radioactive decay of radionuclides and the fourth term describes the mass transfer processes between the components of the solution and formation.

Denote by p the number of chemical reactions involving a j-component, and $${\text{ω}}_i$$ is the reaction rate per unit reactor volume, while the second term on the right side of (2) can be written as

$${\Delta }_j^{ch}={\displaystyle \sum _{i=1}^p{k}_{ji}{\displaystyle \int {\omega }_{i}dV}},$$E3

where $$k_{ji}$$ - the stoichiometric coefficients of chemical reactions of j-component.

The third term can be represented as

$${\Delta }_j^{reac}=-\delta {\displaystyle \int {\left(\alpha \rho \right)}_j}dV,$$E4

where $${\text{α}}_j$$ - the radioactive decay constant of j -th radionuclide, $$\text{δ}$$ - Kronecker delta function, equal to either one if the of j-component of the radionuclide, or zero if otherwise.

We assume that the transition of the impurity molecules of the liquid in the skeleton and its transition from a skeleton into a liquid are determined by the chemical potentials $${\mu }_s,{\mu }_w$$. The fourth term is of the form

$${\Delta }_j^{ex}={\displaystyle \int g\left({\text{μ}}_w{\text{,μ}}_s\right)dV-}{\displaystyle \int g\left({\text{μ}}_j,{\text{μ}}_s\right)dV},$$E5

where $$g\left({\mu }_j,{\mu }_s\right)$$ - a function of mass transfer between the of j-component of the solution and the skeleton of rock, $$g\left({\text{μ}}_w{\text{,μ}}_s\right)$$ - mass transfer function corresponding to the transition of matter from the rock matrix in the solution.

Substituting (3) - (5) into (2) and transforming the surface integral into a volume integral, we obtain

$$\frac{\partial }{\partial \tau }{\displaystyle \int {\rho }_{j}dV}=-{\displaystyle \int \mathrm{div}\left({\overrightarrow{j}}_j\right)dV}+{\displaystyle \sum _{i=1}^p{k}_{ji}{\displaystyle \int {\omega }_{i}dV}}-\delta {\displaystyle \int {\left(\alpha \rho \right)}_j}dV+{\displaystyle \int g\left({\mu }_w,{\mu }_s\right)dV}-{\displaystyle \int g}\left({\mu }_j,{\mu }_s\right)dV.$$E6

By the arbitrariness of the reaction volume $$dV$$ and the continuity of the functions under the integral, we obtain

$$\frac{\partial {\rho }_j}{\partial \tau }=-\mathrm{div}\left({\overrightarrow{j}}_j\right)+{\displaystyle \sum _{i=1}^p{k}_{ji}{\omega }_i}-\delta {\left(\alpha \rho \right)}_j+g\left({\mu }_w,{\mu }_s\right)-g\left({\mu }_j,{\mu }_s\right).$$E7

The resulting equation is nonlinear, even in simple cases the values $${\text{ω}}_i$$ are polynomial functions of concentration. Therefore, in general, equation (7) forms a system of nonlinear partial differential equations. The solution of this system is quite complicated both mathematically, and in terms of its applicability to the description of particular phenomena.

Let us estimate in (7) the contribution of the second term. Obviously, the maximum change of mass in chemical reactions, while other things being equal, will be observed in the following two cases:

$$\begin{array}{l}A+B\to C\uparrow \left(\text{evaporation}\right)\text{,}\\ A+B\to C\downarrow \left(\text{precipitation}\right)\text{.}\end{array}$$

In both reactions the dissolved substances are excluded from consideration, which entails a decrease in the concentrations of the components of the solution. But this type of unpredictable chemical reactions creates the conditions for dangerous situations and in the deep burial of radioactive waste should be excluded. The chemical reaction scheme (acid-base and redox)

$$A+B+\dots \to A^{\prime }+B^{\prime }+\dots \left(\pm \Delta H^0\right).$$E8

are valid and give a slight variation in the concentration of the solution, because typical of enthalpy $$\Delta H^0$$ change is of the order of several hundred kilojoules per mole of interacting substances and the corresponding change in mass $${10}^{-7}÷{10}^{-12}$$ kg, which is negligible in comparison with the mass of dissolved chemical components. Therefore, the change in mass due to chemical reactions will be neglected.

As shown in [35], the time of mass transfer between the fluid and the skeleton of the order of 0.1 s. Thus, the mass transfer, which is characterized by a concentration gradient, is almost instantaneous compared to the time of injection of pollutant that may be from several months to several years. Let us also neglect the processes of chemical compounds leaching from the porous rock to the solution, i.e. assume the condition $$g\left({\text{μ}}_w{\text{,μ}}_s\right)-g\left({\text{μ}}_j,{\text{μ}}_s\right)\approx 0$$.

Based on the above, equation (7) takes the form

$$\frac{\partial {\rho }_j}{\partial \tau }=-\mathrm{div}\left({\overrightarrow{j}}_j\right)-\delta {\left(\alpha \rho \right)}_j.$$E9

Divide the resulting equation into two components: non-radioactive and radioactive fractions

$$\frac{\partial {\rho }_k}{\partial \tau }+\mathrm{div}\left({\overrightarrow{j}}_k\right)=0,\frac{\partial {\rho }_i}{\partial \tau }+\mathrm{div}\left({\overrightarrow{j}}_i\right)=-{\left(\alpha \rho \right)}_i,$$E10

where the index k takes values $$\overline{1,N_{ch}}$$, and the index i takes value $$\overline{1,N_{rc}}$$.

Because of the neglected mass changes in the course of chemical reactions and mass transfer processes in the equilibrium case, it follows that the concentration of impurities non-radioactive fraction with high accuracy can be taken as constant, i. e. $${\rho }_k=\text{const}$$. Then the system of equations (10) can be written as

$$\mathrm{div}\left(\overrightarrow{w}\right)=0,\frac{\partial {\rho }_i}{\partial \tau }+\mathrm{div}\left({\overrightarrow{j}}_i\right)=-{\left(\alpha \rho \right)}_i,$$E11

where $$\text{<mover accent="true">w <mo>→</mo></mover>}$$ - a vector velocity of the fluid.

Write out the flow $${\overrightarrow{j}}_i$$ as the sum of two terms $${\overrightarrow{j}}_i={\overrightarrow{j}}_{Di}+{\rho }_i\overrightarrow{v}$$, where $${\overrightarrow{j}}_{Di}$$ - the diffusion flux, $${\rho }_i\overrightarrow{w}$$- the convective flow, and, taking into account the first equation (11), the second equation takes form

$$\frac{\partial {\rho }_i}{\partial \tau }+\mathrm{div}\left({\overrightarrow{j}}_{Di}\right)+\overrightarrow{w}\cdot \nabla {\rho }_i=-{\left(\alpha \rho \right)}_i.$$E12

According to the Onsager linear theory, the flow for a multicomponent mixture can be written as follows

$${\overrightarrow{j}}_{Di}=-\frac{L_{ij}}{T}{\displaystyle \sum _k\left(\frac{\partial {\mu }_j}{\partial {\rho }_k}\right)}\nabla {\rho }_k,$$E13

where $$L_{ij}$$ - the Onsager kinetic coefficients, $${\text{μ}}_j$$ - the chemical potential of j-th radionuclide.

The real radioactive solutions, arriving at the burial in a deep-seated formations, depending on the half-life have a total volumetric activity of about 10^-6 ~10¹ Ci/l. Let us estimate the mass of radionuclides in solution. Strontium-90 from the volumetric activity 1 Ci/l has a mass of about 7.57 10^-6 kg, and Ruthenium-106 is the same volumetric activity of the mass of the order of 0.3 10^-6 kg. These estimates of the mass of radionuclides provide a basis for considering solution under investigation to be a very dilute solution (with respect to radionuclide fractions). Therefore, the correlation between the diffusion fluxes of components j and k will be negligible.

Thus, the assumption of a very dilute solution leads to the following representation of (13)

$${\overrightarrow{j}}_{Di}=-\frac{L_{ii}}{T}\left(\frac{\partial {\mu }_i}{\partial {\rho }_i}\right)\nabla {\rho }_i,$$E14

Introducing the notation $$D_{ii}=\frac{L_{ii}}{T}\left(\frac{\partial {\mu }_i}{\partial {\rho }_i}\right)$$, we obtain

$${\overrightarrow{j}}_{Di}=-D_{ii}\nabla {\rho }_i.$$E15

Relation (15) is known as Fick's first law, where $$D_{ii}$$ is the diffusion coefficient of i-th - radionuclide.

In many cases, the diffusion coefficient $$D_{ii}$$ can be considered to be constant, then using (15) in equation (12), we obtain a system of equations for evolution of radionuclides in a porous layer

$$\frac{\partial {\rho }_i}{\partial \tau }-D_{ii}\Delta {\rho }_i+\overrightarrow{w}\cdot \nabla {\rho }_i=-{\left(\alpha \rho \right)}_i.$$E16

Equation (16) is written for the porous layer, but it does not take into account the presence of porosity and sorption of radionuclides in the skeleton of the formation. To account for these effects, we introduce an auxiliary space-time function $$m=m\left(t,x,y,z\right)$$, such that $${\displaystyle \int m\left(t,x,y,z\right)dV=V_{por}}$$, where $$V_{por}$$ - the volume of pore space. Obviously, the $$V_s=V-V_{por}$$ - the volume occupied by the formation. Integrating each term of the of equation (16) by volume

$${\displaystyle \int \frac{\partial {\rho }_i}{\partial \tau }}dV-D_{ii}{\displaystyle \int \Delta {\rho }_{i}dV}+\overrightarrow{w}\cdot {\displaystyle \int \nabla {\rho }_{i}dV}=-{\displaystyle \int {\left(\alpha \rho \right)}_i}dV.$$E17

Under the integral expression $${\rho }_{i}dV$$ can be represented as the sum of two terms as $${\rho }_{i}dV={\rho }_{is}d{V}_s+{\rho }_{iw}d{V}_{por}$$, because the other terms, taking into account the mass transfer between the solution and the formation, give a zero contribution due to the steady equilibrium. Therefore

$$\begin{array}{l}{\displaystyle \int \frac{\partial {\rho }_{is}}{\partial \tau }}d{V}_s+{\displaystyle \int \frac{\partial {\rho }_{iw}}{\partial \tau }}d{V}_{por}-D_{is}^i{\displaystyle \int \Delta {\rho }_{is}d{V}_s}-D_{iw}^i{\displaystyle \int \Delta {\rho }_{iw}d{V}_{por}}+\overrightarrow{w}\cdot {\displaystyle \int \nabla {\rho }_{iw}d{V}_{por}}=\\ =-{\displaystyle \int {\left(\alpha {\rho }_s\right)}_i}d{V}_s-{\displaystyle \int {\left(\alpha {\rho }_w\right)}_i}d{V}_{por}.\end{array}$$E18

Using the definition of an auxiliary function $$m$$, it is easy to obtain the following obvious relations

$$d{V}_{por}=mdV,d{V}_s=\left(1-m\right)dV.$$E19

Substituting (19) in equation (18), we obtain

$$\begin{array}{l}{\displaystyle \int (1-m)\frac{\partial {\rho }_{is}}{\partial \tau }}dV+{\displaystyle \int m\frac{\partial {\rho }_{iw}}{\partial \tau }}dV-D_{is}^i{\displaystyle \int (1-m)\Delta {\rho }_{is}dV}-D_{iw}^i{\displaystyle \int m\Delta {\rho }_{iw}dV}+\\ +\overrightarrow{w}\cdot {\displaystyle \int m\nabla {\rho }_{i}dV}=-{\displaystyle \int (1-m){\left(\alpha {\rho }_s\right)}_i}dV-{\displaystyle \int m{\left(\alpha {\rho }_w\right)}_i}dV.\end{array}$$E20

Again, because of the arbitrary choice of the reaction volume $$dV$$ and continuity of integrand functions, we obtain

$$(1-m)\frac{\partial {\rho }_{is}}{\partial \tau }+m\frac{\partial {\rho }_{iw}}{\partial \tau }-D_{is}^i(1-m)\Delta {\rho }_{is}-D_{iw}^{i}m\Delta {\rho }_{iw}+\overrightarrow{w}\cdot m\nabla {\rho }_{iw}=-(1-m){\left(\alpha {\rho }_s\right)}_i-m{\left(\alpha {\rho }_w\right)}_i.$$E21

We assume that the dependence of the impurity concentration in the skeleton of its concentration in the fluid is linear (Henry's isotherm) and does not depend on the volume activity, that is a good approximation for relatively small concentrations of fraction of radionuclide

$${\rho }_{is}=K_{i\Gamma }^i{\rho }_{iw}.$$E22

Then the mass transfer equations take the form:

$$\begin{array}{l}\left[\left(1-m\right)K_{i\Gamma }^i+m\right]\frac{\partial {\rho }_{iw}}{\partial \tau }-\left[D_{is}^i\left(1-m\right)K_{i\Gamma }^i+m{D}_{iw}^i\right]\Delta {\rho }_{iw}+\overrightarrow{w}\cdot m\nabla {\rho }_{iw}=\\ -\left[\left(1-m\right)K_{i\Gamma }^i+m\right]{\left(\alpha {\rho }_w\right)}_i,\end{array}$$E23

where the function $$m$$- void factor, depending on lithological and mineral composition of the layer, $$K_{iГ}^i$$ - the Henry's coefficient of i-th - radionuclide.

The final form of the equations of evolution of radionuclides in solution (liquid phase) in a porous layer, taking into account the porosity and adsorption on the skeleton, one can divide both sides of equation (23) by a factor $$\left(\left(1-m\right)K_{iГ}^i+m\right)$$

$$\frac{\partial {\rho }_{iw}}{\partial \tau }-D_{ii}\Delta {\rho }_{iw}+{\overrightarrow{v^{\prime }}}_{ii}\cdot \nabla {\rho }_{iw}=-{\left(\alpha {\rho }_w\right)}_i,$$E24

here $$D_{ii}=\left(D_{is}^i\left(1-m\right)K_{i\Gamma }^i+m{D}_{iw}^i\right)/\left(\left(1-m\right)K_{i\Gamma }^i+m\right)$$ - the effective diffusion coefficient in the layer, $${{\overrightarrow{v}}^{\prime }}_{ii}=m\overrightarrow{w}/\left(\left(1-m\right)K_{i\Gamma }^i+m\right)$$ - modified velocity of propagation of i-th - radionuclide in a porous layer, (the rate of convective transport of radioactive contaminants).

Note that equation (24) is derived for the case when a radionuclide decaying, forms a non-radioactive nuclide. Possible decay scheme

$$A\to B\to C\left(stable\right),$$

i. e. when the decay product $$B$$ will also be radioactive. The equation takes into account the formation of a child radionuclide

$$\frac{\partial {\rho }_{{}_w}^c}{\partial \tau }-D\Delta {\rho }_{{}_w}^c+\overrightarrow{v^{\prime }}\cdot \nabla {\rho }_{{}_w}^c=\alpha {\rho }_w-{\alpha }^c{\rho }_{{}_w}^c,$$E25

where $${\rho }_{{}_w}^{с}$$ - density of the child radionuclide, $${\rho }_w$$ - density of the parent radionuclide. Investigations of these cases [36], [37] in the work are not included.

The rate of filtration of snap motion of the liquid phases is determined by Darcy's law

$$\overrightarrow{v}=-\frac{k}{\mu }\text{grad}P.$$

In most common filtration processes, the deformation of the porous skeleton, compressibility, and associated with this changes in the temperature of liquids are small. The main effects that determine the motion of the system are the non-equlibrium joint motion of several liquid phases, molecular and convective diffusion of solute in the phases of the components, the absorption of the solid phase or sorption of the components, mass transfer between phases.

Thus, the system of equations describing the mass transfer during injection of liquid radioactive wastes in deep porous horizon is as follows:

$$\begin{array}{l}\mathrm{div}\left(\overrightarrow{w}\right)=0,\\ \frac{\partial {\rho }_{iw}}{\partial \tau }-D_i\Delta {\rho }_{iw}+{\overrightarrow{v^{\prime }}}_i\cdot \nabla {\rho }_{iw}=-{\alpha }_i{\rho }_{iw},\\ \frac{\partial {\rho }_{{}_w}^{с}}{\partial \tau }-D\Delta {\rho }_{{}_w}^{с}+\overrightarrow{v^{\prime }}\cdot \nabla {\rho }_{{}_w}^{с}=\alpha {\rho }_w-{\alpha }^{с}{\rho }_{{}_w}^{с}.\end{array}$$E26

For a complete statement of the problem requires knowledge of the radiochemical composition of the solution, flow rate of the injection, diffusion parameters and the geometry of the simulated porous layer. Note that if the injection rate is known, it is easy to determine the rate of filtration. Then integrating Darcy's equation, we can describe the pressure field in the formation.

The problem under consideration has cylindrical symmetry about the axis of the well, through which the liquid wastes are ejected; it is convenient to represent the system of equations (26) in a cylindrical coordinate system.

Writing the first equation (26) in a cylindrical coordinate system and, given that the liquid is distributed in the porous layer only in the radial direction, we obtain the equation for the velocity field:

$$\frac{\partial \left(r{\text{w}}_r\right)}{\partial r_{\text{d}}}=0,$$

solving this equation and applying the obvious boundary condition $${{\text{w}}_r|}_{r=r_0}={\text{w}}_{\text{0}}$$, where $${\text{w}}_{\text{0}}$$ - velocity of the fluid from the cased hole in the porous layer, we have

$${\text{w}}_r={\text{w}}_0{r}_0/r_{\text{d}}.$$E27

Then the remaining equations of (26) using (27), and the anisotropy of diffusion coefficients and thermal conductivity in the directions r and z in a cylindrical coordinate system can be written as

$$\begin{array}{l}\frac{\partial {\rho }_{iw}}{\partial \tau }+\frac{{v^{\prime }}_0{}_i{r}_0}{r_d}\frac{\partial {\rho }_{iw}}{\partial r_d}-D_r{}_i\frac{1}{r_d}\frac{\partial }{\partial r_d}\left(r_d\frac{\partial {\rho }_{iw}}{\partial r_d}\right)-D_z{}_i\frac{{\partial }^2{\rho }_{iw}}{\partial z_d{}^2}=-{\alpha }_i{\rho }_{iw},\\ \frac{\partial {\rho }_w^{с}}{\partial \tau }+\frac{{v^{\prime }}_0{r}_0}{r_d}\frac{\partial {\rho }_w^{с}}{\partial r_F}-D_r\frac{1}{r_d}\frac{\partial }{\partial r_d}\left(r_d\frac{\partial {\rho }_w^c}{\partial r_d}\right)-D_z\frac{{\partial }^2{\rho }_w^c}{\partial z_d{}^2}=\alpha {\rho }_w-\alpha {\rho }_w^c,\end{array}$$E28

where $${\text{v}}_{\text{0}}=m{\text{w}}_{\text{0}}$$ - the rate of fluid filtration, $$D_{ri},D_{zi}$$ - an effective diffusion coefficient in the direction r_d and z_d, respectively, and the multiplier $${\text{γ}}_i=\left[\left(1-m\right)K_{i\Gamma }+m\right]$$.

It is assumed that the real porous layer is represented by a multiphase system, where each phase consists of a sufficiently large number of randomly distributed small particles. Particle size, small in comparison with the basic physical quantities are assumed to be so large that within each particle condition of "local equilibrium" [38] and all the conservation laws are satisfied [9]. All contact surfaces of particles of different nature are surfaces of discontinuity of some physical fields. However, the above assumptions allow us in physically small volumes to define the space of continuous functions, carrying out the description of the fields of each phase. This determination is carried out by a predetermined method of averaging, from which, in general, depend on the results obtained in [9]. As with most occurring filtration processes, the deformation of the porous skeleton, compressibility and associated changes in temperature fluids rely small.

Given that the determining factor in the process of mass transfer is the concentration of the parent nuclide, confine ourselves to the problem for a single pollutant, which is radioactive and chemically active. The first equation (28) is represented as

$$\frac{\partial {\rho }_w}{\partial \tau }+\frac{\overrightarrow{v}\nabla {\rho }_w}{\left(1-m\right)K_{\Gamma }+m}-D\Delta {\rho }_w=-\alpha {\rho }_w.$$E29

Here we have introduced the notation

$$D=\frac{D_s\left(1-m\right)K_{\text{г}}+D_{w}m}{\left(1-m\right)K_{\text{г}}+m}$$E30

D is the effective diffusion coefficient in the layer. From (29) that in the equation describing the migration of contaminants, it is necessary to take into account the convective transport of pollutants, "complicated by" the presence of porosity in the skeleton and mass transfer processes occurring between the pollutant and the skeleton. Equation (29) to determine the rate of convective transport of pollutants in porous layer $${\overrightarrow{v}}^{\prime }$$, by analogy with the rate of convective heat transfer and flow rate $$\overrightarrow{v}$$

$${\overrightarrow{v}}^{\prime }=\frac{\overrightarrow{v}}{\left(1-m\right)K_{\Gamma }+m}.$$E31

The rate of convective transport of the impurity $${\overrightarrow{v}}^{\prime }$$ determines the position of the front of pollution R_d, just as the filtration rate $$\overrightarrow{v}$$ determines the position of the front of injected fluid R_w. The position of the injected fluid front is determined from the mass balance of the injected fluid and, for the case of injection at a constant speed v₀ into the layer through a cased hole of radius r₀, the corresponding expression is given by

$$R_w=\sqrt{\frac{2v_0{r}_0\tau }{m}+r_0^2}=\sqrt{2w_0{r}_0\tau +r_0^2}=\sqrt{\frac{Q\tau }{\pi mH}+r_0^2}.$$E32

3. The mathematical formulation of the problem of mass transfer

Fig. 1 shows the geometry of the problem in a cylindrical coordinate system whose axis coincides with the axis of the borehole. The environment is presented by three areas with flat boundaries. Injection of impurities into the area is out of the hole radius, covering and underlying layers are impermeable, middle area is a the porous region, all layers are considered homogeneous and anisotropic on the diffusion properties. Observation is carried out at a distance from the axis of the borehole

Figure 1.

The geometry of the problem: 0, 1, 2 - porous, covering and underlying layers, respectively, 3 – borehole

Through a hole of small (compared to the distance to the observation point) radius $$r_0$$ in an infinite horizontal layer of thickness $$-h<z_{\text{d}}<h$$ water with a radioactive contaminant are injected. In arriving liquid at $$r\le r_0$$ the concentration of impurities kept constant and equal to $${\rho }_0$$. The concentration of pollutants in the layer changes due to convective transport along the direction $$r$$, the diffusion along $$r,z$$ and the concentration of sources. As such sources of radioactive decay of pollutant are considered. Field of densities during the filtration of radioactive solutions was investigated in [36, 39-50].

The mathematical formulation of the problem of mass transfer for all areas involves the diffusion equation with taking into account the radioactive decay in the covering

$$\frac{\partial {\rho }_{1d}}{\partial \tau }-D_{\text{z1}}\frac{{\partial }^2{\rho }_{1d}}{\partial z_d^2}-D_{\text{r1}}\frac{1}{r_d}\frac{\partial }{\partial r_d}\left(r_d\frac{\partial {\rho }_{1d}}{\partial r_d}\right)=-\text{α}{\rho }_{1d},\tau >0,r_d>0,z_d>h$$E33

and the underlying

$$\frac{\partial {\rho }_{2d}}{\partial \tau }-D_{z2}\frac{{\partial }^2{\rho }_{2d}}{\partial z_d^2}-D_{r2}\frac{1}{r_d}\frac{\partial }{\partial r_d}\left(r_d\frac{\partial {\rho }_{2d}}{\partial r_d}\right)=-\text{α}{\rho }_{2d},\tau >0,r_d>0,z_d<-h$$E34

layers, as well as the equation of convective diffusion, taking into account the radioactive decay in the porous layer

$$\frac{\partial {\rho }_d}{\partial \tau }-D_z\frac{{\partial }^2{\rho }_d}{\partial z_d^2}-D_r\frac{1}{r_d}\frac{\partial }{\partial r_d}\left(r_d\frac{\partial {\rho }_d}{\partial r_d}\right)+\frac{{\text{<msup>v <mo>′</mo></msup>}}_{\text{0}}{r}_{\text{0}}}{r_{\text{d}}}\frac{\partial {\rho }_d}{\partial r_d}=-\text{α}{\rho }_d,\tau >0,r_d>0,\left|z_d\right|<h.$$E35

The conditions of conjugation represent the equality of densities and fluxes of dissolved substances at the interface of the layers

$${{\rho }_d|}_{z_d=h}={{\rho }_{d1}|}_{z_d=h},{{\rho }_d|}_{z_d=-h}={{\rho }_{d2}|}_{z_d=-h},$$E36

$${D_z\frac{\partial {\rho }_d}{\partial z_d}|}_{z_d=h}={D_{z1}\frac{\partial {\rho }_{1d}}{\partial z_d}|}_{z_d=h},{D_z\frac{\partial {\rho }_d}{\partial z_d}|}_{z_d=-h}={D_{z2}\frac{\partial {\rho }_{2d}}{\partial z_d}|}_{z_d=-h}.$$E37

The density of pollutant at the entrance of porous layer assumed to be constant

$${{\rho }_d|}_{r_d=0}=\left[m+K\left(1-m\right)\right]{\rho }_0.$$E38

Assuming that at the initial time the density of the of pollutant is equal to zero

$${{\rho }_d|}_{\tau =0}={{\rho }_{1d}|}_{\tau =0}={{\rho }_{2d}|}_{\tau =0}=0.$$E39

In addition, at infinity the conditions of regularity

$${{\rho }_d|}_{r_d\to +\infty }=0,{{\rho }_{1d}|}_{r_d+z_d\to +\infty }=0,{{\rho }_{2d}|}_{r_d+\left|z_d\right|\to +\infty }=0.$$E40

Let us turn then to the dimensionless quantities

$$\begin{array}{l}\rho =\frac{{\rho }_d}{{\rho }_0},r=\frac{r_d}{h},z=\frac{z_d}{h},t=\frac{{\lambda }_{z1}\tau }{c_1{\rho }_{п1}{h}^2},\mathrm{At}=\frac{c_1{\rho }_{п1}}{{\lambda }_{z1}}\alpha h^2,{}_1^{2}D=\frac{D_{z2}}{D_{z1}},\\ {}_1^{0}D=\frac{D_z}{D_{z1}},{}_2^{0}D=\frac{D_z}{D_{z2}},a_{z1}=\frac{{\lambda }_{z1}}{c_1{\rho }_{п1}},\gamma =\frac{1}{m+K\left(1-m\right)}.\end{array}$$

We also introduce the analogue of the Péclet number

$$\mathrm{Pd}={v^{\prime }}_0{r}_0/D_{z1},$$

where $${v^{\prime }}_0$$ is the rate of convective transport of the pollutant at a distance $$r_0$$ from the axis of the borehole. With this notation equations (33) - (40) take the form

$$\frac{\partial {\rho }_1}{\partial t}-\frac{{\partial }^2{\rho }_1}{\partial z^2}-\frac{D_{r1}}{D_{z1}}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial {\rho }_1}{\partial r}\right)=-\text{At}{\rho }_1,t>0,r>0,z>1,$$E41

$$\frac{\partial {\rho }_{\text{2}}}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2}{\partial z^2}-\frac{D_{r2}}{D_{z1}}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial {\rho }_2}{\partial r}\right)=-\text{At}{\rho }_2,t>0,r>0,z<-1,$$E42

$$\frac{\partial \rho }{\partial t}-{}_1^{0}D\frac{{\partial }^2\rho }{\partial z^2}-\frac{D_r}{D_{z1}}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}\right)+\frac{\text{Pd}}{r}\frac{\partial \rho }{\partial r}=-\text{At}\rho ,t>0,r>0,\left|z\right|<1.$$E43

Let us estimate the ratio of the third and fourth terms in equation (43)

$$\frac{\frac{D_r}{D_{z1}}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}\right)}{\frac{\text{Pd}}{r}\frac{\partial \rho }{\partial r}}\approx \frac{D_r\frac{{\rho }_0}{R^2}}{D_{z1}\text{Pd}\frac{{\rho }_0}{R^2}}=\frac{D_r}{D_{z1}\text{Pd}}\approx \frac{1}{\text{Pd}}<<1,$$

Boundary, initial conditions and conjugation conditions are not changed

$${\rho |}_{t=0}={{\rho }_1|}_{t=0}={{\rho }_2|}_{t=0}=0,$$E44

$${\rho |}_{r=0}=1/\gamma ,{\rho |}_{z=1}={{\rho }_1|}_{z=1},{\rho |}_{z=-1}={{\rho }_2|}_{z=-1},$$E45

$${\frac{\partial \rho }{\partial z}|}_{z=1}={{}_0^{1}D\frac{\partial {\rho }_1}{\partial z}|}_{z=1},{\frac{\partial \rho }{\partial z}|}_{z=-1}={{}_0^{2}D\frac{\partial {\rho }_2}{\partial z}|}_{z=-1},$$E46

$${\rho |}_{r\to \infty }=0,{{\rho }_1|}_{r+z\to \infty }=0,{{\rho }_2|}_{r+\left|z\right|\to \infty }=0.$$E47

The system of equations (44) - (50) defines a mathematical formulation of the problem of mass transfer.

4. Expansion of the solution to the problem of mass transfer on the asymptotic parameter

Let us consider the more general problem, which is obtained by introducing into the equations and boundary conditions of arbitrary asymptotic parameter $$\epsilon$$ of formal substitution in the diffusion coefficient $$D_z$$ for $$D_z/\text{ε}$$. In accordance with the designations this performed by replacing $${}_0^{1}D$$ for $$\text{ε}{}_0^{1}D$$ and $${}_0^{2}D={}_1{}^{2}D{}_0^{1}D$$ for $$\text{ε}{}_0^{2}D$$. Note that the original problem can be obtained from the solution of a parameterized the problem when $$\epsilon =1$$. The problem (44) - (50) is thus a particular case of the more general parameterized problem, containing a parameter of the asymptotic expansion both in the equation for the layer and in the conditions of conjugation

$$\frac{\partial {\rho }_1}{\partial t}-\frac{{\partial }^2{\rho }_1}{\partial z^2}=-\text{At}{\rho }_1,t>0,r>0,z>1,$$E51

$$\epsilon \frac{\partial \rho }{\partial t}{-}_1^{0}D\frac{{\partial }^2\rho }{\partial z^2}+\epsilon \frac{\mathrm{Pd}}{r}\frac{\partial \rho }{\partial r}=-\epsilon \text{At}\rho ,t>0,r>0,\left|z\right|<1,$$E52

$$\frac{\partial {\rho }_2}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2}{\partial z^2}=-\text{At}{\rho }_2,t>0,r>0,z<-1$$E53

with boundary conditions

$${\frac{\partial \rho }{\partial z}|}_{z=1}={{}_0^{1}D\epsilon \frac{\partial {\rho }_1}{\partial z}|}_{z=1},{\frac{\partial \rho }{\partial z}|}_{z=-1}={}_0^{2}D{\epsilon \frac{\partial {\rho }_2}{\partial z}|}_{z=-1},$$E54

$${\rho |}_{r=0}=1/\gamma ,{\rho |}_{z=1}={{\rho }_1|}_{z=1},{\rho |}_{z=-1}={{\rho }_2|}_{z=-1},$$E55

$${\rho |}_{t=0}={{\rho }_1|}_{t=0}={{\rho }_2|}_{t=0}=0,$$E56

$${\rho |}_{r\to +\infty }=0,{{\rho }_1|}_{r+z\to +\infty }=0,{{\rho }_2|}_{r+\left|z\right|\to +\infty }=0.$$E57

To find the solution to (51) - (57), one can represent the density function $$\rho$$ of each region by the asymptotic formula of the parameter $$\epsilon$$

$$\rho ={\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\mathrm{...}+{\epsilon }^n{\rho }^{\left(n\right)}+{\theta }^{\left(n\right)},{\rho }_i={\rho }_i^{(0)}+\epsilon {\rho }_i^{(1)}+\mathrm{...}+{\epsilon }^n{\rho }_i^{\left(n\right)}+{\theta }_i^{\left(n\right)},i=1,2.$$E58

Substituting expression (58) in (51) - (57) and grouping terms in powers of the expansion parameter $$\text{ε}$$, one can easily obtain

$$\left(\frac{\partial {\rho }_1^{(0)}}{\partial t}-\frac{{\partial }^2{\rho }_1^{(0)}}{\partial z^2}+\mathrm{At}{\rho }_1^{\left(0\right)}\right)+\epsilon \left(\frac{\partial {\rho }_1^{(1)}}{\partial t}-\frac{{\partial }^2{\rho }_1^{(1)}}{\partial z^2}+\mathrm{At}{\rho }_1^{\left(1\right)}\right)+\mathrm{...}=0,t>0,r>0,z>1,$$E59

$$\left(\frac{\partial {\rho }_2^{\left(0\right)}}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2^{\left(0\right)}}{\partial z^2}+\mathrm{At}{\rho }_2^{\left(0\right)}\right)+\epsilon \left(\frac{\partial {\rho }_2^{\left(1\right)}}{\partial t}{-}_1^{2}D\frac{{\partial }^2{\rho }_2^{\left(1\right)}}{\partial z^2}+\mathrm{At}{\rho }_2^{\left(1\right)}\right)+\mathrm{...}=0,t>0,r>0,z<-1,$$E60

$$\begin{array}{l}{-}_1^{0}D\frac{{\partial }^2{\rho }^{\left(0\right)}}{\partial z^2}+\epsilon \left(\frac{\partial {\rho }^{\left(0\right)}}{\partial t}-{}_1^{0}D\frac{{\partial }^2{\rho }^{\left(1\right)}}{\partial z^2}+\frac{Pd}{r}\frac{\partial {\rho }^{\left(0\right)}}{\partial r}+\mathrm{At}{\rho }^{\left(0\right)}\right)+\\ +{\epsilon }^2\left(\frac{\partial {\rho }^{\left(1\right)}}{\partial t}-{}_1^{0}D\frac{{\partial }^2{\rho }^{\left(2\right)}}{\partial z^2}+\frac{Pd}{r}\frac{\partial {\rho }^{\left(1\right)}}{\partial r}+\mathrm{At}{\rho }^{\left(1\right)}\right)+\mathrm{...}=0,t>0,r>0,\left|z\right|<1,\end{array}$$E61

$${\frac{\partial {\rho }^{\left(0\right)}}{\partial z}|}_{z=1}+\epsilon \left({\frac{\partial {\rho }^{\left(1\right)}}{\partial z}|}_{z=1}-{{}_0^{1}D\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}\right)+\mathrm{...}=0,{\frac{\partial {\rho }^{\left(0\right)}}{\partial z}|}_{z=-1}+\epsilon \left({\frac{\partial {\rho }^{\left(1\right)}}{\partial z}|}_{z=-1}-{{}_0^{2}D\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right)+\mathrm{...}=0,$$E62

$${\left({\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\mathrm{...}\right)|}_{z=1}={\left({\rho }_1^{\left(0\right)}+\epsilon {\rho }_1^{\left(1\right)}+\dots \right)|}_{z=1},{\left({\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\mathrm{...}\right)|}_{z=-1}=\left({\rho }_2^{\left(0\right)}+\epsilon {\rho }_2^{\left(1\right)}+\mathrm{...}\right)|{}_{z=-1},$$E63

$${\left({\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\dots \right)|}_{t=0}={\left({\rho }_1^{\left(0\right)}+\epsilon {\rho }_1^{\left(1\right)}+\dots \right)|}_{t=0}={\left({\rho }_2^{\left(0\right)}+\epsilon {\rho }_2^{\left(1\right)}+\dots \right)|}_{t=0}=0,$$E64

$${\left({\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\mathrm{...}\right)|}_{r=0}=1/\gamma ,$$E65

$$\left({\rho }^{\left(0\right)}+\epsilon {\rho }^{\left(1\right)}+\dots \right)|{}_{r\to +\infty }=0,{\left({\rho }_1^{\left(0\right)}+\epsilon {\rho }_1^{\left(1\right)}+\dots \right)|}_{r+z\to +\infty }=0,{\left({\rho }_2^{\left(0\right)}+\epsilon {\rho }_2^{\left(1\right)}+\dots \right)|}_{r+\left|z\right|\to +\infty }=0.$$E66

Analysis of the formulation of the problem shows that the factors of powers of ε in (61) contain the neighboring coefficients of the expansion, and in this sense, are linked. To solve the corresponding equations implemented decoupling procedure.

4.1. The mathematical formulation of the problem of mass transfer in a zero approximation

If we formally consider ε in equation (61) to be infinitely small, we obtain $${}_1^{0}D{\partial }^2{\rho }^{\left(0\right)}/\partial z^2=0$$. The result of integration $$\partial {\rho }^{(0)}/\partial z=A\left(r,t\right)$$ with the boundary conditions (62) allows us to establish that $$A\left(r,t\right)=0.$$Thus, in the zero approximation, the density of the pollutant $${\rho }^{\left(0\right)}={\rho }^{\left(0\right)}\left(r,t\right)$$ is a function only of r and t. Consequently, in the zero approximation the density of the pollutant in each cylindrical cross section with the axis $$z$$ is the same in height of the carrier layer. Next, equating to zero coefficients near $$\epsilon$$ in equation (61), we obtain

$$\frac{\partial {\rho }^{(0)}}{\partial t}-{}_1^{0}D\frac{{\partial }^2{\rho }^{\left(1\right)}}{\partial z^2}+\frac{Pd}{r}\frac{\partial {\rho }^{\left(0\right)}}{\partial r}+\mathrm{At}{\rho }^{\left(0\right)}=0.$$E67

Since $${\rho }^{\left(0\right)}\left(r,t\right)$$ does not depend on $$z,$$ the auxiliary function $$E\left(r,t\right)$$, composed of the terms of the equation (67) containing $${\text{ρ}}^{\left(0\right)}$$

$$E\left(r,t\right)=\frac{\partial {\rho }^{(0)}}{\partial t}+\frac{Pd}{r}\frac{\partial {\rho }^{\left(0\right)}}{\partial r}+\mathrm{At}{\rho }^{\left(0\right)},$$E68

is also independent of z. Then (1) can be written as

$${}_1^{0}D\frac{{\partial }^2{\rho }^{\left(1\right)}}{\partial z^2}=E\left(r,t\right).$$E69

Integrating successively, one can find the expression for the first derivative of the first coefficient $${\rho }^{\left(1\right)}$$ of the variable $$\text{z}$$

$$\frac{\partial {\rho }^{\left(1\right)}}{\partial z}{=}_0^{1}D\left[zE\left(r,t\right)+F\left(r,t\right)\right],$$E70

and the first coefficient of expansion in the form of a quadratic trinomial

$${\rho }^{\left(1\right)}={}_0^{1}D\left(\frac{z^{\text{2}}}{2}E\left(r,t\right)+zF\left(r,t\right)+Q\left(r,t\right)\right),$$E71

with the functional coefficients to be determined. From the boundary conditions (62) with the cofactor $$\epsilon$$ we have

$${{}_1^{2}D\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}=E\left(r,t\right)+F\left(r,t\right),{{}_1^{2}D\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}=-E\left(r,t\right)+F\left(r,t\right).$$E72

Hence, one can obtain an expression for the functional coefficients and through the traces of derivatives of the outer regions

$$E\left(r,t\right)=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right),$$E73

$$F\left(r,t\right)=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}+{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right).$$E74

Substituting (73) in (74), one can obtain the desired equation for the zero approximation of the density of impurities in the layer

$$\frac{\partial {\rho }^{\left(0\right)}}{\partial t}+\frac{\mathrm{Pd}}{r}\frac{\partial {\rho }^{\left(0\right)}}{\partial t}+\mathrm{At}{\rho }^{\left(0\right)}=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right).$$E75

The final statement of the problem in the zero approximation also includes the equations in the covering and underlying formations

$$\frac{\partial {\rho }_1^{\text{(0)}}}{\partial t}-\frac{{\partial }^2{\rho }_1^{\text{(0)}}}{\partial z^2}=-\text{At}{\rho }_1^{\left(0\right)},t>0,r>0,z>1,$$E76

$$\frac{\partial {\rho }_2^{\text{(0)}}}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2^{\text{(0)}}}{\partial z^2}=-\text{At}{\rho }_2^{\left(0\right)},t>0,r>0,z<-1,$$E77

$$\frac{\partial {\rho }^{\left(0\right)}}{\partial t}+\frac{\mathrm{Pd}}{r}\frac{\partial {\rho }^{\left(0\right)}}{\partial r}+\mathrm{At}{\rho }^{\left(0\right)}=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right),t>0,r>0,\left|z\right|<1,$$E78

and the appropriate initial and boundary conditions

$${{\rho }^{\left(0\right)}|}_{t=0}={{\rho }_1^{(0)}|}_{t=0}={{\rho }_2^{(0)}|}_{t=0}=0,$$E79

$${\rho }^{\left(0\right)}={{\rho }_1^{(0)}|}_{z=1}={{\rho }_2^{\left(0\right)}|}_{z=-1},$$E80

$${{\rho }^{\left(0\right)}|}_{r=0}=1/\gamma ,$$E81

$${{\rho }^{\left(0\right)}|}_{r\to +\infty }=0,{{\rho }_1^{(0)}|}_{r+z\to +\infty }=0,{{\rho }_2^{(0)}|}_{r+\left|z\right|\to +\infty }=0.$$E82

Expressions (76) - (82) represent the boundary value problem for zero expansion coefficient $${\rho }^{\left(0\right)}$$ or zero approximation. Note that in contrast to the original, which is the problem of conjugation for parabolic equations, it is mixed, since the equation contains traces of derivatives from the outer regions.

Finding the zero approximation of the density of the radioactive contaminant is important because just that approach arises in the zero approximation of the temperature problem.

4.2. The zero approximation of the problem of mass transport as the solution of the averaged problem

Let us average a parameterized problem (51) - (57) over z within the carrier layer according to

$$⟨\rho ⟩=\frac{1}{2}{\displaystyle \underset{-1}{\overset{1}{\int }}\rho dz}.$$

Successively averaging each term of equation (34)

$$\begin{array}{l}⟨\frac{\partial \rho }{\partial t}⟩=\frac{1}{2}{\displaystyle \underset{-1}{\overset{1}{\int }}\frac{\partial \rho }{\partial t}dz}=\frac{\partial ⟨\rho ⟩}{\partial t},\mathrm{Pd}\frac{1}{r}⟨\frac{\partial \rho }{\partial r}⟩=\mathrm{Pd}\frac{1}{r}\frac{\partial ⟨\rho ⟩}{\partial r}\\ \frac{1}{\epsilon }\cdot {}_1^{0}D⟨\frac{{\partial }^2\rho }{\partial z^2}⟩=\frac{{}_1^{0}D}{2\epsilon }{\displaystyle \underset{-1}{\overset{1}{\int }}\frac{{\partial }^2\rho }{\partial z^2}dz=}\frac{{}_1^{0}D}{2\epsilon }\left({\frac{\partial \rho }{\partial z}|}_{z=1}-{\frac{\partial \rho }{\partial z}|}_{z=-1}\right)=\frac{1}{2}\left({\frac{\partial {\rho }_1}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2}{\partial z}|}_{z=-1}\right),\end{array}$$

one can obtain the following formulation of the averaged parameterized problem:

$$\frac{\partial {\rho }_1}{\partial t}-\frac{{\partial }^2{\rho }_1}{\partial z^2}+\mathrm{At}{\rho }_1=0,t>0,r>0,z>1,$$E83

$$\frac{\partial ⟨\rho ⟩}{\partial t}-\frac{1}{2}\left({\frac{\partial {\rho }_1}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2}{\partial z}|}_{z=-1}\right)+\mathrm{Pd}\frac{1}{r}\frac{\partial ⟨\rho ⟩}{\partial r}=-\mathrm{At}⟨\rho ⟩,t>0,r>0,\left|z\right|<1,$$E84

$$\frac{\partial {\rho }_2}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2}{\partial z^2}+\mathrm{At}{\rho }_2=0,t>0,r>0,z<-1$$E85

$$⟨\rho ⟩={{\rho }_1|}_{z=1}={{\rho }_2|}_{z=-1},$$E86

$${⟨\rho ⟩|}_{r=0}=1/\gamma ,$$E87

$${⟨\rho ⟩|}_{t=0}={{\rho }_1|}_{t=0}={{\rho }_2|}_{t=0}=0,$$E88

$${⟨\rho ⟩|}_{r\to +\infty }=0,{{\rho }_1|}_{r+z\to +\infty }=0,{{\rho }_2|}_{r+\left|z\right|\to +\infty }=0.$$E89

The resulting problem coincides with problem (83) - (89) for the zero approximation of the density of pollutant. The uniqueness of solutions implies that $$⟨\rho ⟩={\rho }^{\left(0\right)}$$, i. e. zero approximation describes the known way averaged solutions to the original problem. If average the original nonparameterized problem (44) - (50), it also coincides with the problem for zero approximation of the field densities of the pollutant.

4.3. The mathematical formulation of the problem of mass transfer in the first approximation

Equations (41) - (43) for the coefficients for $$\epsilon$$ take the form

$$\frac{\partial {\rho }_1^{(1)}}{\partial t}-\frac{{\partial }^2{\rho }_1^{(1)}}{\partial z^2}+\mathrm{At}{\rho }_1^{\left(1\right)}=0,t>0,r>0,z>1,$$E90

$$\frac{\partial {\rho }_2^{(1)}}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\rho }_2^{(1)}}{\partial z^2}+\mathrm{At}{\rho }_2^{\left(1\right)}=0,t>0,r>0,z<-1,$$E91

$$\frac{\partial {\rho }^{\left(1\right)}}{\partial t}+\frac{\mathrm{Pd}}{r}\frac{\partial {\rho }^{\left(1\right)}}{\partial r}-{}_1^{0}D\frac{{\partial }^2{\rho }^{\left(2\right)}}{\partial z^2}+\mathrm{At}{\rho }^{\left(1\right)}=0,t>0,r>0,\left|z\right|<1,$$E92

appropriate boundary and initial conditions are represented as

$${\frac{\partial {\rho }^{\left(1\right)}}{\partial z}|}_{z=1}-{{}_0^{1}D\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}=0,{\frac{\partial {\rho }^{\left(1\right)}}{\partial z}|}_{z=-1}-{{}_0^{2}D\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}=0,$$E93

$${{\rho }^{\left(1\right)}|}_{z=1}={{\rho }_1^{(1)}|}_{z=1},{{\rho }^{\left(1\right)}|}_{z=-1}={{\rho }_2^{\left(1\right)}|}_{z=-1},$$E94

$${{\rho }^{\left(1\right)}|}_{t=0}={{\rho }_1^{\left(1\right)}|}_{t=0}={{\rho }_2^{\left(1\right)}|}_{t=0}=0,$$E95

$${{\rho }^{\left(1\right)}|}_{r\to +\infty }=0,{{\rho }_1^{\left(1\right)}|}_{r+z\to +\infty }=0,{{\rho }_2^{\left(1\right)}|}_{r+\left|z\right|\to +\infty }=0,$$E96

$${{\rho }^{\left(1\right)}|}_{r=0}=0.$$E97

The solution to the problem is sought in the form of quadratic polynomial for z according to (71), where the auxiliary functions $$E\left(r,t\right)$$ and $$F\left(r,t\right)$$ are defined by (73) - (74), and the function $$Q\left(r,t\right)$$ is unknown. For its determination we can write (92) as

$$\frac{{\partial }^2{\rho }^{\left(2\right)}}{\partial z^2}=\frac{{}_0^{1}D}{\delta }\stackrel{⌢}{L}{\rho }^{\left(1\right)},$$E98

here $$\stackrel{⌣}{L}$$ is the operator

$$\stackrel{⌢}{L}=\frac{\partial }{\partial t}+\frac{\text{Pd}}{r}\frac{\partial }{\partial r}+\text{At}.$$E99

Using (71) and (98), and linearity of the operator $$\stackrel{⌣}{L}$$, one can obtain

$$\frac{{\partial }^2{\rho }^{\left(2\right)}}{\partial z^2}=\frac{{}_0^1{D}^2}{{\delta }^2}\left(\frac{z^2}{2}\stackrel{⌢}{L}E\left(r,t\right)+z\stackrel{⌢}{L}F\left(r,t\right)+\stackrel{⌢}{L}Q\left(r,t\right)\right).$$E100

Integrating the last expression over z

$$\frac{\partial {\rho }^{\left(2\right)}}{\partial z}=\frac{{}_0^1{D}^2}{{\delta }^2}\left(\frac{z^3}{6}\stackrel{⌢}{L}E\left(r,t\right)+\frac{z^2}{2}\stackrel{⌢}{L}F\left(r,t\right)+z\stackrel{⌢}{L}Q\left(r,t\right)+W\left(r,t\right)\right).$$E101

From the expression (101) and the boundary conditions (62) we have

$${{}_1^{0}D\frac{\partial {\rho }_1^{\left(1\right)}}{\partial z}|}_{z=1}=\frac{1}{6}\stackrel{⌢}{L}E\left(r,t\right)+\frac{1}{2}\stackrel{⌢}{L}F\left(r,t\right)+\stackrel{⌢}{L}Q\left(r,t\right)+W\left(r,t\right),$$E102

$${}_1^{0}D{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)}}{\partial z}|}_{z=-1}=-\frac{1}{6}\stackrel{⌢}{L}E\left(r,t\right)+\frac{1}{2}\stackrel{⌢}{L}F\left(r,t\right)-\stackrel{⌢}{L}Q\left(r,t\right)+W\left(r,t\right).$$E103

From (90) and (91) one can get the equation for the definition of Q

$$\stackrel{⌢}{L}Q\left(r,t\right)=\frac{{}_1^{0}D}{2}\left({\frac{\partial {\rho }_1^{\left(1\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)}}{\partial z}|}_{z=-1}\right)-\frac{1}{6}\stackrel{⌢}{L}E\left(r,t\right).$$E104

The equation for the determination of first coefficient of expansion is obtained by substituting (57), (58) and (55) in (98) with using (100)

$$\begin{array}{l}\stackrel{⌢}{L}{\rho }^{\left(1\right)}=\frac{{}_0^{1}D}{4}\left(z^2-\frac{1}{3}\right)\stackrel{⌢}{L}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right)+\\ \frac{{}_0^{1}D}{2}z\stackrel{⌢}{L}\left({\frac{\partial {\rho }_1^{\left(0\right)}}{\partial z}|}_{z=1}+{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)}}{\partial z}|}_{z=-1}\right)+\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(1\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)}}{\partial z}|}_{z=-1}\right).\end{array}$$E105

The problem to determine the first coefficient of expansion also includes equations (90), (91) for the density field in the covering and underlying layers, respectively.

It is easy to see that the problem formulated by (90), (91), (105), (93) - (96) with the condition (97) has a trivial solution, so the condition (97) is redundant and should be weakened.

4.4. The problem for the remainder — Additional boundary condition

Substituting the expansion (58) with n = 1 in the parameterized problem (51) - (57), we obtain a problem for the remainder term

$$\frac{\partial {\theta }_1}{\partial t}-\frac{{\partial }^2{\theta }_1}{\partial z^2}=-\text{At}{\theta }_1,t>0,r>0,z>1,$$E106

$$\epsilon \frac{\partial \theta }{\partial t}{-}_1^{0}D\frac{{\partial }^2\theta }{\partial z^2}+\epsilon \frac{\mathrm{Pd}}{r}\frac{\partial \theta }{\partial r}+\epsilon \mathrm{At}\theta =-{\epsilon }^2\stackrel{⌢}{L}{\rho }^{\left(1\right)},t>0,r>0,\left|z\right|<1,$$E107

$$\frac{\partial {\theta }_2}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\theta }_{\text{2}}}{\partial z^2}=-\text{At}{\theta }_2,t>0,r>0,z<-1$$E108

with the boundary conditions and conjugation conditions

$${\frac{\partial \theta }{\partial z}|}_{z=1}={}_0^{1}D{\epsilon \frac{\partial {\theta }_1}{\partial z}|}_{z=1}{+{}_0^{1}D{\epsilon }^2\frac{\partial {\rho }_1^{(1)}}{\partial z}|}_{z=1},{\frac{\partial \theta }{\partial z}|}_{z=-1}={}_0^{2}D{\epsilon \frac{\partial {\theta }_2}{\partial z}|}_{z=-1}{+{\text{}}_0^{2}D{\epsilon }^2\frac{\partial {\rho }_2^{(1)}}{\partial z}|}_{z=-1},$$E109

$${\theta |}_{z=1}={{\theta }_1|}_{z=1},{\theta |}_{z=-1}={{\theta }_2|}_{z=-1},$$E110

$${\theta |}_{t=0}={{\theta }_1|}_{t=0}={{\theta }_2|}_{t=0}=0,$$E111

$${\theta |}_{r=0}=-{\epsilon {\rho }^{\left(1\right)}|}_{r=0},$$E112

$${\theta |}_{r\to +\infty }=0,{{\theta }_1|}_{r+z\to +\infty }=0,{{\theta }_2|}_{r+\left|z\right|\to +\infty }=0.$$E113

Restrict our investigation of the problem, averaged over the thickness of the layer. By averaging the second derivative over the vertical coordinate, use the conjugation conditions (109)

$$\begin{array}{l}\frac{{\text{}}_1^{0}D}{\epsilon }⟨\frac{{\partial }^2\theta }{\partial z^2}⟩=\frac{{\text{}}_1^{0}D}{2\epsilon }{\displaystyle \underset{-1}{\overset{1}{\int }}\frac{{\partial }^2\theta }{\partial z^2}dz=}\frac{{\text{}}_1^{0}D\delta }{2\epsilon }\left({\frac{\partial \theta }{\partial z}|}_{z=1}-{\frac{\partial \theta }{\partial z}|}_{z=-1}\right)=\\ \frac{1}{2}\left({\frac{\partial {\theta }_1}{\partial z}|}_{z=1}-{{}_1^{2}D\frac{\partial {\theta }_2}{\partial z}|}_{z=-1}\right)+\frac{\epsilon }{2}\left({\frac{\partial {\rho }_1^{(1)}}{\partial z}|}_{z=1}{-{\text{}}_1^{2}D\frac{\partial {\rho }_2^{(1)}}{\partial z}|}_{z=-1}\right).\end{array}$$

The final formulation of the averaged problem for the remainder term represented as

$$\frac{\partial {\theta }_1}{\partial t}-\frac{{\partial }^2{\theta }_1}{\partial z^2}=-\text{At}{\theta }_1,t>0,r>0,z>1,$$E114

$$\begin{array}{l}\frac{\partial ⟨\theta ⟩}{\partial t}+\frac{Pd}{r}\frac{\partial ⟨\theta ⟩}{\partial r}+\text{At}⟨\theta ⟩-\frac{1}{2}\left({\frac{\partial {\theta }_1}{\partial z}|}_{z=1}-{{}_1^{2}D\frac{\partial {\theta }_2}{\partial z}|}_{z=-1}\right)=\\ =-\text{ε}\left(\stackrel{⌢}{L}⟨{\rho }^{\left(1\right)}⟩-\frac{1}{2}\left({\frac{\partial {\rho }_1^{(1)}}{\partial z}|}_{z=1}{-{\text{}}_1^{2}D\frac{\partial {\rho }_2^{(1)}}{\partial z}|}_{z=-1}\right)\right),t>0,r>0,\left|z\right|<1,\end{array}$$E115

$$\frac{\partial {\theta }_2}{\partial t}-{}_1^{2}D\frac{{\partial }^2{\theta }_{\text{2}}}{\partial z^2}=-\text{At}{\theta }_2,t>0,r>0,z<-1,$$E116

$${⟨\theta ⟩={\theta }_1|}_{z=1}={{\theta }_2|}_{z=-1},$$E117

$${⟨\theta ⟩|}_{r=0}=-{\epsilon ⟨{\rho }^{\left(1\right)}⟩|}_{r=0},$$E118

$${⟨\theta ⟩|}_{t=0}={{\theta }_1|}_{t=0}={{\theta }_2|}_{t=0}=0,$$E119

$$⟨\theta ⟩|{}_{r\to +\infty }=0,{{\theta }_1|}_{r+z\to +\infty }=0,{{\theta }_2|}_{r+\left|z\right|\to +\infty }=0$$E120

It is easy to show that the averaged problem (110) - (116) for the remainder term has a trivial solution if and only if

$$⟨{\rho }^{\left(1\right)}⟩|{}_{r=0}=0,$$E121

$$\stackrel{⌢}{L}⟨{\rho }^{\left(1\right)}⟩-\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(1\right)}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)}}{\partial z}|}_{z=-1}\right)=0,$$E122

that is, when in the averaged problem for remainder term there are no sources. Averaging (107) with regard to (109), one can show that (122) is satisfied identically. Thus, in order to the averaged problem for the remainder term had a trivial solution it is sufficient for the average condition (121) to be satisfied. Therefore, in order to get the exact on the average solution of (44) - (50) on the field of density in the layer, in the formulation of the problem for the first coefficients of the asymptotic expansion (90) (91), (105) (95) - (97) the boundary condition (97) must be replaced by non-local (121).

5. Solution to the problem of mass transfer in the zero approximation

5.1. Solution to the problem in the zero approximation

In the image space of Laplace-Carson the problem (76) - (82) in the zero approximation is represented as

$$p{\rho }_1^{\left(0\right)u}-\frac{{\partial }^2{\rho }_1^{\left(0\right)u}}{\partial z^2}=-\mathrm{At}{\rho }_1^{\left(0\right)u},r>0,z>1,$$E123

$$p{\rho }_2^{\left(0\right)u}-{}_1^{2}D\frac{{\partial }^2{\rho }_2^{\left(0\right)u}}{\partial z^2}=-\mathrm{At}{\rho }_2^{\left(0\right)u},r>0,z<-1,$$E124

$$p{\rho }^{\left(0\right)u}+\frac{\text{Pd}}{r}\frac{\partial {\rho }^{\left(0\right)u}}{\partial r}+\text{At}{\rho }^{\left(0\right)u}=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}\right)\text{,}r>0\text{,}\left|z\right|<1\text{,}$$E125

$${{\rho }_1^{\left(0\right)u}|}_{z=1}={\rho }^{\left(0\right)u}={{\rho }_2^{\left(0\right)u}|}_{z=-1},$$E126

$${{\rho }^{\left(0\right)u}|}_{r=0}=1/\gamma ,$$E127

$${{\rho }^{\left(0\right)u}|}_{r\to \infty }=0,{{\rho }_1^{\left(0\right)u}|}_{r+z\to \infty }=0,{{\rho }_2^{\left(0\right)u}|}_{r+\left|z\right|\to \infty }=0.$$E128

Taking into account the boundary conditions (126), as well as the fact that in the zero approximation the density of the pollutant in the porous layer is independent of z and is a function only of r and t, the solution of equations (123), (124) can be rewritten as follows:

$${\rho }_1^{\left(0\right)u}={\rho }^{\left(0\right)u}\exp \left(-\sqrt{p+\mathrm{At}}\left(z-1\right)\right),$$E129

$${\rho }_2^{\left(0\right)u}={\rho }^{\left(0\right)u}\exp \left(\sqrt{{}_2^{1}D\left(p+\mathrm{At}\right)}\left(z+1\right)\right).$$E130

These expressions allow us to determine the values of the traces of derivatives from the outer regions included in the equation for the layer, through the density of impurities in it

$${\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}=-\sqrt{p+\mathrm{At}}{\rho }^{\left(0\right)u},{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}=\sqrt{{}_2^{1}D\left(p+\mathrm{At}\right)}{\rho }^{\left(0\right)u}.$$E131

Substituting (131) into equation (125), after simple transformations we obtain an ordinary differential equation for the determination of $${\rho }^{\left(0\right)u}$$

$$\frac{\mathrm{Pd}}{r}\frac{d{\rho }^{\left(0\right)u}}{dr}=-\left(p+\mathrm{At}+\frac{1}{2}\sqrt{p+\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)\right){\rho }^{\left(0\right)u},$$E132

from which we finally get

$${\rho }^{\left(0\right)u}=\gamma \exp \left[-\left(p+\mathrm{At}+\frac{1}{2}\sqrt{p+\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)\right)\frac{r^2}{2\mathrm{Pd}}\right].$$E133

Translation in the original space is carried out by the reference [51]. The expression for the density of radioactive contaminants for the porous layer in the original space is represented as

$$\begin{array}{l}{\rho }^{\left(0\right)}=\frac{\gamma }{2}\text{exp}\left(-\text{At}{r}^2/2\text{Pd}\right)\Phi \left(t-r^2/2\text{Pd}\right)\times \\ \times [\text{exp}\left(\frac{-\sqrt{\text{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2}{4\text{Pd}}\right)\text{erfc}\left(\frac{\sqrt{\delta }\left(1+\sqrt{{}_1^{2}D}\right)r^2}{8\text{Pd}\sqrt{t-r^2/2\text{Pd}}}-\sqrt{\text{At}\left(t-r^2/2\text{Pd}\right)}\right)+\\ +\text{exp}\left(\frac{\sqrt{\text{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2}{4\text{Pd}}\right)\text{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2}{8\text{Pd}\sqrt{t-r^2/2\text{Pd}}}+\sqrt{\text{At}\left(t-r^2/2\text{Pd}\right)}\right)],t>0,r>0,\left|z\right|<1.\end{array}$$E134

Also transition is feasible in the original space for the coating (129) and underlying (130) layers

$$\begin{array}{l}{\rho }_1^{\left(0\right)}=\frac{\gamma }{2}\exp \left(-\mathrm{At}{r}^2/2\mathrm{Pd}\right)\Phi \left(t-r^2/2\mathrm{Pd}\right)\times \\ \times \{\exp \left(-\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2/4\mathrm{Pd}-\sqrt{\mathrm{At}}\left(z-1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2+4\mathrm{Pd}\left(z-1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}-\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)+\\ +\exp \left(\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2/4\mathrm{Pd}+\sqrt{\mathrm{At}}\left(z-1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2+4\mathrm{Pd}\left(z-1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}+\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)\},t>0,r>0,z>1,\end{array}$$E135

$$\begin{array}{l}{\rho }_2^{\left(0\right)}=\frac{\gamma }{2}\exp \left(-\mathrm{At}{r}^2/2\mathrm{Pd}\right)\Phi \left(t-r^2/2\mathrm{Pd}\right)\times \{\exp \left(-\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2/4\mathrm{Pd}+\sqrt{\mathrm{At}{}_2^{1}D}\left(z+1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2-4\mathrm{Pd}\sqrt{{}_2^{1}D}\left(z+1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}-\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)+\\ +\exp \left(\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)r^2/4\mathrm{Pd}-\sqrt{{\mathrm{At}}_2^{1}D}\left(z+1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2-4\mathrm{Pd}\sqrt{{}_2^{1}D}\left(z+1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\gamma Pt}}+\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)\},t>0,r>0,z<-1.\end{array}$$E136

The first factor in the solution (134) - (136) describes the decrease in the density of the pollutant as a result of radioactive decay, the second - the Heaviside function, determines the radius of the spread of contamination zone and the third (the expression in curly brackets) takes into account changes of the density due to the diffusion of pollutants, and radioactive decay of diffusing nuclide. Since the contribution of radioactive decay is described by the factor of $$\text{exp}\left(-\text{At}{r}^2/2\text{Pd}\right)$$, then it can be argued that the concentration of the radioactive contaminant is reduced by a factor of e due to the decay at the distances defined by the simple relation R_e=h$$\sqrt{2\text{Pd/At}}$$=$$\sqrt{{\text{2v}}_{\text{0}}{r}_0/\alpha }$$. It follows that for short-lived isotopes, zone of contamination is low. On the other hand, to reduce the zone of influence of long-lived radioactive isotopes, the rate of filtration should be reduced.

The resulting solution (134) contains the Heaviside function, which vanishes for r ≥$$\sqrt{2\text{Pd}t},$$and helps to determine the radius of the zone of radioactive contamination

$$R_p=h\sqrt{2\text{Pd}t=}\sqrt{{\text{2v}}_{\text{0}}{r}_0\tau }.$$E137

The most important physical results are described by the zero approximation of the asymptotic expansion, the first and the following coefficients determine the "correction term". In addition, due to the smallness of the diffusion coefficient ($$D_z$$~10^-9÷10^-11), the spread of a contaminant in water-resistant layers in the vertical direction is negligible compared to the convective transport in a porous layer and has little effect on the size of the zone of contamination.

5.2. The solution to the problem of mass transfer in the first approximation

In the space transformations of Laplace-Carson, the problem (90), (91), (105), (95) - (97), (121) for the first coefficient of expansion is represented as

$$\begin{array}{l}{\stackrel{⌢}{L}}^u{\rho }^{\left(1\right)u}=\frac{{}_0^{1}D}{4}\left(z^2-\frac{1}{3}\right){\stackrel{⌢}{L}}^u\left({\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}\right)+\\ +\frac{{}_0^{1}D}{2}z{\stackrel{⌢}{L}}^u\left({\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}+{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}\right)+\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(1\right)u}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)u}}{\partial z}|}_{z=-1}\right),r>0,\left|z\right|<1,\end{array}$$E138

$$p{\rho }_1^{\left(1\right)u}-\frac{{\partial }^2{\rho }_1^{\left(1\right)u}}{\partial z^2}+\mathrm{At}{\rho }_1^{\left(1\right)u}=0,r>0,z>1,$$E139

$$p{\rho }_2^{\left(1\right)u}-{}_1^{2}D\frac{{\partial }^2{\rho }_2^{\left(1\right)u}}{\partial z^2}+\mathrm{At}{\rho }_2^{\left(1\right)u}=0,r>0,z<-1.$$E140

The initial conditions and conjugation conditions at the boundaries are represented as

$${{\rho }^{\left(1\right)u}|}_{z=1}={{\rho }_1^{(1)u}|}_{z=1},$$E141

$${{\rho }^{\left(1\right)u}|}_{z=-1}={{\rho }_2^{\left(1\right)u}|}_{z=-1},$$E142

$${{\rho }^{\left(1\right)u}|}_{r\to +\infty }=0,{{\rho }_1^{\left(1\right)u}|}_{r+z\to +\infty }=0,{{\rho }_2^{\left(1\right)u}|}_{r+\left|z\right|\to +\infty }=0,$$E143

$${⟨{\rho }^{\left(1\right)u}⟩|}_{r=0}=0.$$E144

The operator $$\stackrel{⌣}{L}$$ in the image space has the form

$${\stackrel{⌣}{L}}^u=p+\text{At}+\gamma \frac{\text{Pt}}{r}\frac{\partial }{\partial r}.$$

The action of this operator on the zero expansion in the image space is determined by the formula

$${\stackrel{⌢}{L}}^u{\rho }^{\left(0\right)u}=-\frac{\sqrt{p+\text{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right){\rho }^{\left(0\right)u}.$$E145

The solution in the first approximation, according to (71), is sought in the form of quadratic polynomial

$${\rho }^{\left(1\right)u}={}_0^{1}D\left(\frac{z^2}{2}{E}^u+z{F}^u+Q^u\right),$$E146

in which $$E^u$$ and $$F^u$$ are expressed through the zero approximation according to (73), (74)

$$E^u=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}\right),F^u=\frac{1}{2}\left({\frac{\partial {\rho }_1^{\left(0\right)u}}{\partial z}|}_{z=1}+{}_1^{2}D{\frac{\partial {\rho }_2^{\left(0\right)u}}{\partial z}|}_{z=-1}\right),$$E147

and the function $$Q^u$$, according to (104), is defined by the equation

$${\stackrel{⌢}{L}}^u{Q}^u=\frac{{}_1^{0}D}{2}\left({\frac{\partial {\rho }_1^{\left(1\right)u}}{\partial z}|}_{z=1}-{}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)u}}{\partial z}|}_{z=-1}\right)-\frac{1}{6}{\stackrel{⌢}{L}}^u{E}^u.$$E148

Solutions of equations (139), (140) shall be as follows:

$${\rho }_1^{\left(1\right)u}={{\rho }^{\left(1\right)u}|}_{z=1}\exp \left(-\sqrt{p+\mathrm{At}}\left(z-1\right)\right),{\rho }_2^{\left(1\right)u}={{\rho }^{\left(1\right)u}|}_{z=-1}\exp \left(\sqrt{{}_2^{1}D\left(p+\mathrm{At}\right)}\left(z+1\right)\right).$$E149

Let us find the traces of the outer regions of the right side of equation (138)

$${\frac{\partial {\rho }_1^{\left(1\right)u}}{\partial z}|}_{z=1}=-{}_0^{1}D\sqrt{p+\mathrm{At}}\left(\frac{1}{2}{E}^u+F^u+Q^u\right),$$E150

$${}_1^{2}D{\frac{\partial {\rho }_2^{\left(1\right)u}}{\partial z}|}_{z=-1}{=}_0^{2}D\sqrt{{}_2^{1}D\left(p+\mathrm{At}\right)}\left(\frac{1}{2}{E}^u-F^u+Q^u\right).$$E151

Note also that the action of the operator $${\stackrel{⌣}{L}}^{и}$$ on functions

$$E^u=-\frac{\sqrt{p+\mathrm{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right){\rho }^{\left(0\right)u},F^u=-\frac{\sqrt{\delta \left(p+\mathrm{At}\right)}}{2}\left(1-\sqrt{{}_1^{2}D}\right){\rho }^{\left(0\right)u}$$E152

according to (147), this leads to the following:

$${\stackrel{⌢}{L}}^u{E}^u=\frac{p+\mathrm{At}}{4}{\left(1+\sqrt{{}_1^{2}D}\right)}^2{\rho }^{\left(0\right)u},{\stackrel{⌢}{L}}^u{F}^u=\frac{p+\mathrm{At}}{4}{\left(1-{}_1^{2}D\right)}^2{\rho }^{\left(0\right)u}.$$E153

The final equation for determination of $$Q^{и}$$ takes the form

$$\frac{\mathrm{Pd}}{r}\frac{d{Q}^u}{dr}+\left(p+\mathrm{At}+\frac{\sqrt{p+\mathrm{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right)\right)Q^u=\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{3}\left(p+\mathrm{At}\right){\rho }^{\left(0\right)u}.$$E154

It’s general solution is represented as

$$\begin{array}{l}{Q}^u=\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{3\mathrm{Pd}}\delta \left(p+\mathrm{At}\right)\times {\displaystyle \underset{0}{\overset{r}{\int }}{\rho }^{\left(0\right)u}\exp \left[-\left(p+\mathrm{At}+\frac{\sqrt{p+\mathrm{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right)\right)\frac{r^2-{r^{\prime }}^2}{2\mathrm{Pd}}\right]r^{\prime }d{r}^{\prime }+}\\ \text{}\text{}+C\exp \left[-\left(p+\mathrm{At}+\frac{\sqrt{p+\mathrm{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right)\right)\frac{r^2}{2\mathrm{Pd}}\right].\end{array}$$E155

Constant C is determined from the averaging (144)

$${Q^u|}_{r=0}=-\frac{1}{6}{E^u|}_{r=0}.$$E156

Hence, one can obtain

$$Q^u=\frac{1}{6}\left[\frac{\left(1+{}_1^{2}D-\sqrt{{}_1^{2}D}\right)}{\mathrm{Pd}}\delta \left(p+\mathrm{At}\right)r^2+\frac{\sqrt{p+\mathrm{At}}}{2}\left(1+\sqrt{{}_1^{2}D}\right)\right]{\rho }^{\left(0\right)u}.$$E157

As a result, the solution for the first coefficient of expansion in the images is represented as

$${\rho }^{\left(1\right)u}=\frac{{}_0^{1}D}{2}\sqrt{p+\mathrm{At}}[\left(\frac{1}{6}-\frac{z^2}{2}\right)\left(1+\sqrt{{}_1^{2}D}\right)-z\left(1-\sqrt{{}_1^{2}D}\right)]{\rho }^{\left(0\right)u}+{}_0^{1}D\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{6\mathrm{Pd}}{r}^2\left(p+\mathrm{At}\right){\rho }^{\left(0\right)u},$$E158

$${\rho }_1^{\left(1\right)u}=\frac{{}_0^{1}D}{3}\left[\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{2Pd}{r}^2\left(p+\mathrm{At}\right)-\sqrt{p+\mathrm{At}}\left(2-\sqrt{{}_1^{2}D}\right)\right]{\rho }^{\left(0\right)u}\times \exp \left(-\left(z-1\right)\sqrt{p+\mathrm{At}}\right),$$E159

$${\rho }_2^{\left(1\right)u}=\frac{{}_0^{1}D}{3}\left[\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{2\mathrm{Pd}}{r}^2\left(p+\mathrm{At}\right)+\sqrt{p+\mathrm{At}}\left(1-2\sqrt{{}_1^{2}D}\right)\right]{\rho }^{\left(0\right)u}\times \exp \left(\left(z+1\right)\sqrt{{}_2^{1}D\left(p+\mathrm{At}\right)}\right).$$E160

Determination of the originals is carried out by help of the following correspondence:

$$\begin{array}{l}\sqrt{p+\gamma }\exp \left(-\sqrt{\beta \left(p+\gamma \right)}\right)\to \frac{1}{\sqrt{\pi t}}\exp \left(-\gamma t-\frac{\beta }{4t}\right)-\\ -\frac{\sqrt{\gamma }}{2}[\exp \left(\sqrt{\beta \gamma }\right)\mathrm{erfc}\left(\frac{1}{2}\sqrt{\frac{\beta }{t}}+\sqrt{\gamma t}\right)-\exp \left(-\sqrt{\beta \gamma }\right)\mathrm{erfc}\left(\frac{1}{2}\sqrt{\frac{\beta }{t}}-\sqrt{\gamma t}\right)],\\ p\exp \left(-\sqrt{\beta \left(p+\gamma \right)}\right)\to \frac{\sqrt{\beta }}{2t\sqrt{\pi t}}\exp \left(-\frac{\beta }{4t}-\gamma t\right).\end{array}$$

Finally, we obtain for the porous layer

$$\begin{array}{l}{\rho }^{\left(1\right)}=\frac{{}_0^{1}D}{4\gamma }\sqrt{\mathrm{At}}\exp \left(-\mathrm{At}{r}^2/2\mathrm{Pd}\right)\Phi \left(t-r^2/2\mathrm{Pd}\right)\times \\ \times [\frac{2\left(1+\sqrt{{}_1^{2}D}\right)\sqrt{\mathrm{At}\left(t-r^2/2\mathrm{Pd}\right)}}{\sqrt{\pi }}(\frac{1-\sqrt{{}_1^{2}D}+{}_1^{2}D}{24{\mathrm{Pd}}^2\left(t-r^2/2\mathrm{Pd}\right)}{r}^4+\left(\frac{1}{6}-\frac{z^2}{2}-z\frac{1-\sqrt{{}_1^{2}D}}{1+\sqrt{{}_1^{2}D}}\right))\times \\ \times \exp \left(-\frac{{\left(1+\sqrt{{}_1^{2}D}\right)}^2}{64{\mathrm{Pd}}^2\left(t-r^2/2\mathrm{Pd}\right)}{r}^4-\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)\right)+\\ +\left(\sqrt{\mathrm{At}}\frac{1-\sqrt{{}_1^{2}D}+{}_1^{2}D}{3\mathrm{Pd}}{r}^2-\left(1+\sqrt{{}_1^{2}D}\right)\left(\frac{1}{6}-\frac{z^2}{2}-z\frac{1-\sqrt{{}_1^{2}D}}{1+\sqrt{{}_1^{2}D}}\right)\right)\times \\ \times \exp \left(\frac{\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\mathrm{Pd}}{r}^2\right)\mathrm{erfc}\left(\frac{1+\sqrt{{}_1^{2}D}}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}{r}^2+\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)+\\ +\left(\sqrt{\mathrm{At}}\frac{1-\sqrt{{}_1^{2}D}{+}_1^{2}D}{3\mathrm{Pd}}{r}^2+\left(1+\sqrt{{}_1^{2}D}\right)\left(\frac{1}{6}-\frac{z^2}{2}-z\frac{1-\sqrt{{}_1^{2}D}}{1+\sqrt{{}_1^{2}D}}\right)\right)\times \\ \times \exp \left(-\frac{\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\mathrm{Pd}}{r}^2\right)\mathrm{erfc}\left(\frac{1+\sqrt{{}_1^{2}D}}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}{r}^2-\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)],r>0,t>0,\left|z\right|<1,\end{array}$$E161

for the covering layer

$$\begin{array}{l}{\rho }_1^{\left(1\right)}={}_0^{1}D\frac{\sqrt{\text{At}}}{6\gamma }\text{exp}\left(-\frac{\text{At}}{\text{2Pd}}{r}^2\right)\Phi \left(t-\frac{r^2}{2\text{Pd}}\right)\times \\ \times [\frac{2}{\sqrt{\pi }\sqrt{\text{At}\left(t-r^{\text{2}}/2\text{Pd}\right)}}\left(\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{4\text{Pd}\left(t-r^{\text{2}}/2\text{Pd}\right)}{r}^2\left(\frac{1+\sqrt{{}_1^{2}D}}{4\text{Pd}}{r}^2+\text{z}-1\right)-\left(2-\sqrt{{}_1^{2}D}\right)\right)\times \\ \times \text{exp}\left(-\frac{{\left(\left(1+\sqrt{{}_1^{2}D}\right)r^2+4\text{Pd}\left(\text{z}-\text{1}\right)\right)}^2}{64{\text{Pd}}^{\text{2}}\left(t-r^{\text{2}}/2\text{Pd}\right)}-\text{At}\left(t-\frac{r^{\text{2}}}{2\text{Pd}}\right)\right)+\\ +\left(\frac{\sqrt{\text{At}}\left(1+{}_1^{2}D-\sqrt{{}_1^{2}D}\right)}{2\text{Pd}}{r}^2-\left(2-\sqrt{{}_1^{2}D}\right)\right)\text{exp}\left(-\frac{\sqrt{\text{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\text{Pd}}{r}^2-\sqrt{\text{At}}\left(\text{z}-1\right)\right)\times \\ \times \text{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2}{8\text{Pd}\sqrt{t-r^{\text{2}}/2\text{Pd}}}+\frac{\text{z}-\text{1}}{2\sqrt{t-r^{\text{2}}/2\text{Pd}}}-\sqrt{\text{At}\left(t-\frac{r^{\text{2}}}{2\text{Pd}}\right)}\right)+\\ +\left(\frac{\sqrt{At}\left(1+{}_1{}^{2}D-\sqrt{{}_1{}^{2}D}\right)}{2\text{Pd}}{r}^2+\left(2-\sqrt{{}_1{}^{2}D}\right)\right)\text{exp}\left(-\frac{\sqrt{\text{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\text{Pd}}{r}^2+\sqrt{\text{At}}\left(\text{z}-1\right)\right)\times \\ \times \text{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2}{8\text{Pd}\sqrt{t-r^{\text{2}}/2\text{Pd}}}+\frac{\text{z}-\text{1}}{2\sqrt{t-r^{\text{2}}/2\text{Pd}}}+\sqrt{\text{At}\left(t-\frac{r^{\text{2}}}{2\text{Pd}}\right)}\right).\end{array}$$E162

and for the underlying layer

$$\begin{array}{l}{\rho }_2^{\left(1\right)}={}_0^{1}D\frac{\sqrt{\mathrm{At}}}{6\gamma }\exp \left(-\frac{\mathrm{At}}{2\mathrm{Pd}}{r}^2\right)\Phi \left(t-\frac{r^2}{2\mathrm{Pd}}\right)\times \times [\frac{2}{\sqrt{\pi }\sqrt{\mathrm{At}\left(t-r^2/2\mathrm{Pd}\right)}}\left(\frac{1+{}_1^{2}D-\sqrt{{}_1^{2}D}}{4\mathrm{Pd}\left(t-r^2/2\mathrm{Pd}\right)}{r}^2\left(\frac{1+\sqrt{{}_1^{2}D}}{4\mathrm{Pd}}{r}^2-\sqrt{{}_2^{1}D}\left(z+1\right)\right)+1-2\sqrt{{}_1^{2}D}\right)\times \\ \times \exp \left(-\frac{{\left(\left(1+\sqrt{{}_1^{2}D}\right)r^2-4\mathrm{Pd}\sqrt{{}_2^{1}D}\left(z+1\right)\right)}^2}{64{\mathrm{Pd}}^2\left(t-r^2/2\mathrm{Pd}\right)}-\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)\right)+\\ +\left(\frac{\sqrt{\mathrm{At}}\left(1+{}_1^{2}D-\sqrt{{}_1^{2}D}\right)}{2\mathrm{Pd}}{r}^2+1-2\sqrt{{}_1^{2}D}\right)\exp \left(-\frac{\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\mathrm{Pd}}{r}^2+\sqrt{{}_2^{1}D\mathrm{At}}\left(z+1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2-4\mathrm{Pd}\sqrt{{}_2^{1}D}\left(z+1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}-\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)+\\ +\left(\frac{\sqrt{\mathrm{At}}\left(1+{}_1^{2}D-\sqrt{{}_1^{2}D}\right)}{2\mathrm{Pd}}{r}^2-1+2\sqrt{{}_1^{2}D}\right)\exp \left(\frac{\sqrt{\mathrm{At}}\left(1+\sqrt{{}_1^{2}D}\right)}{4\mathrm{Pd}}{r}^2-\sqrt{{}_2^{1}D\mathrm{At}}\left(z+1\right)\right)\times \\ \times \mathrm{erfc}\left(\frac{\left(1+\sqrt{{}_1^{2}D}\right)r^2-4\mathrm{Pd}\sqrt{{}_2^{1}D}\left(z+1\right)}{8\mathrm{Pd}\sqrt{t-r^2/2\mathrm{Pd}}}+\sqrt{\mathrm{At}\left(t-\frac{r^2}{2\mathrm{Pd}}\right)}\right)],r>0,t>0,z<-1.\end{array}$$E163

Note that when r = 0 the first coefficient of expansion (158)

$${\rho }^{\left(1\right)u}|{}_{r=0}=\frac{{}_0^{1}D}{2}\sqrt{p+\mathrm{At}}[\frac{1}{2}\left(\frac{1}{3}-z^2\right)\left(1+\sqrt{{}_1^{2}D}\right)-z\left(1-\sqrt{{}_1^{2}D}\right)]$$E164

depends on the variable z and does not vanish, that is constructed solution does not satisfy the boundary condition (48). To eliminate the marked lack, it is necessary to build a boundary-layer functions in a neighborhood of r = 0.

5.3. Construction of boundary-layer solutions

Let us represent the solution to (78) - (82) as

$${\rho }_1={\stackrel{⌢}{\rho }}_1+{\Pi }_1,\rho =\stackrel{⌢}{\rho }+\Pi ,{\rho }_2={\stackrel{⌢}{\rho }}_2+{\Pi }_2$$E165

where $$\stackrel{⌢}{\rho }=\stackrel{⌢}{\rho }\left(r,z,t\right)$$ - the regular part, $$\Pi =\Pi \left(y,z,t\right)$$ - boundary-layer part of the expansion by the asymptotic parameter, $$y=r^2/2\text{ε}$$ the stretched variable [44]. Substituting (165) to (44) - (50) and applying the transformation of Laplace-Carson over the variable t, in the image space we obtain a problem for the boundary-layer functions

$${\Pi }^{и}={\Pi }^{\left(0\right)u}+\text{ε}{\Pi }^{\left(1\right)u}+\mathrm{...},{\Pi }_1^u={\Pi }_1^{\left(0\right)u}+\text{ε}{\Pi }_1^{\left(1\right)u}+\mathrm{...},{\Pi }_2^u={\Pi }_2^{\left(0\right)u}+\text{ε}{\Pi }_2^{\left(1\right)u}+\mathrm{...}$$E166

in the following way:

$$\left(\frac{{\partial }^2{\Pi }_1^{\left(0\right)u}}{\partial z^2}-\left(p+\mathrm{At}\right){\Pi }_1^{\left(0\right)u}\right)+\epsilon \left(\frac{{\partial }^2{\Pi }_1^{\left(1\right)u}}{\partial z^2}-\left(p+\mathrm{At}\right){\Pi }_1^{\left(1\right)u}\right)+\mathrm{...}=0,y>0,z>1,$$E167

$$\left(\frac{{\partial }^2{\Pi }_{}^{\left(0\right)u}}{\partial z^2}{-}_0^{1}D\mathrm{Pd}\frac{\partial {\Pi }_{}^{\left(0\right)u}}{\partial y}\right)+\epsilon \left(\frac{{\partial }^2{\Pi }_{}^{\left(1\right)u}}{\partial z^2}{-}_0^{1}D\mathrm{Pd}\frac{\partial {\Pi }_{}^{\left(1\right)u}}{\partial y}-\left(p+\mathrm{At}\right){\Pi }_{}^{\left(0\right)u}\right)+\mathrm{...}=0,y>0,\left|z\right|<1,$$E168

$$\left(\frac{{\partial }^2{\Pi }_2^{\left(0\right)u}}{\partial z^2}{-}_2^{1}D\left(p+\mathrm{At}\right){\Pi }_2^{\left(0\right)u}\right)+\epsilon \left(\frac{{\partial }^2{\Pi }_2^{\left(1\right)u}}{\partial z^2}{-}_2^{1}D\left(p+\mathrm{At}\right){\Pi }_2^{\left(1\right)u}\right)+\mathrm{...}=0,y>0,z<-1,$$E169

$$\frac{\partial {\Pi }^{\left(0\right)u}}{\partial z}|{}_{z=1}+\epsilon \left(\frac{\partial {\Pi }^{\left(1\right)u}}{\partial z}|{}_{z=1}-{}_0^{1}D\frac{\partial {\Pi }_1^{\left(0\right)u}}{\partial z}|{}_{z=1}\right)+\mathrm{...}=0,\frac{\partial {\Pi }^{\left(0\right)u}}{\partial z}|{}_{z=-1}+\epsilon \left(\frac{\partial {\Pi }^{\left(1\right)u}}{\partial z}|{}_{z=-1}-{}_0^{2}D\frac{\partial {\Pi }_2^{\left(0\right)u}}{\partial z}|{}_{z=-1}\right)+\mathrm{...}=0,$$E170

$$\left({\Pi }^{\left(0\right)и}-{\Pi }_1^{\left(0\right)u}+\epsilon \left({\Pi }^{\left(1\right)u}-{\Pi }_1^{\left(1\right)u}\right)+\mathrm{...}\right)|{}_{z=1}=0,\left({\Pi }^{\left(0\right)u}-{\Pi }_2^{\left(0\right)u}+\epsilon \left({\Pi }^{\left(1\right)u}-{\Pi }_2^{\left(1\right)u}\right)+\mathrm{...}\right)|{}_{z=-1}=0,$$E171

$${\Pi }^{\left(0\right)u}|{}_{y=0}+\epsilon \left({\rho }^{\left(1\right)u}+{\Pi }^{\left(1\right)u}\right)|{}_{y=0}+\mathrm{...}=0,$$E172

$$\left({\Pi }^{\left(0\right)u}+\epsilon {\Pi }^{\left(1\right)u}+\mathrm{...}\right)|{}_{y\to \infty }=0,\left({\Pi }_1^{\left(0\right)u}+\epsilon {\Pi }_1^{\left(1\right)u}+\mathrm{...}\right)|{}_{y+z\to \infty }=0,\left({\Pi }_2^{\left(0\right)u}+\epsilon {\Pi }_2^{\left(1\right)u}+\mathrm{..}\right)|{}_{y+\left|z\right|\to \infty }=0.$$E173

The problem for the zero coefficients has only the trivial solution

$${\Pi }^{\left(0\right)u}={\Pi }_1^{\left(0\right)u}={\Pi }_2^{\left(0\right)u}=0.$$E174

The problem for the first coefficient is divided into three independent parts. The first one is the problem for the boundary layer functions in the layer of the form

$$\frac{{\partial }^2{\Pi }^{\left(1\right)u}}{\partial z^2}-\left(p+\mathrm{At}\right)\frac{\partial {\Pi }^{\left(1\right)u}}{\partial y}=0,y>0,\left|z\right|<1,$$E175

$$\frac{\partial {\Pi }^{\left(1\right)u}}{\partial z}|{}_{z=1}=\frac{\partial {\Pi }^{\left(1\right)u}}{\partial z}|{}_{z=-1}=0,$$E176

$${\Pi }^{\left(1\right)u}|{}_{y=0}=-{\rho }^{\left(1\right)u}|{}_{y=0}=-\sqrt{p+\mathrm{At}}\frac{{}_0^{1}D}{2}\left[\left(\frac{1}{6}-\frac{z^2}{2}\right)\left(1+\sqrt{{}_1^{2}D}\right)-z\left(1-\sqrt{{}_1^{2}D}\right)\right],$$E177

$${\Pi }^{\left(1\right)u}|{}_{y\to \infty }=0.$$E178

The solution to this problem found by the method of separation of variables and is defined by the formula

$$\begin{array}{l}{\Pi }^{\left(1\right)u}=\frac{4\sqrt{p+\mathrm{At}}{}_1^{0}D}{{\pi }^2}\left(1-\sqrt{{}_1^{2}D}\right)\sin \left(\frac{\pi }{2}z\right)\exp \left(-{\left(\frac{\pi }{2}\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right)+\\ +\frac{\sqrt{p+\mathrm{At}}{}_1^{0}D}{{\pi }^2}{\displaystyle \sum _{n=1}^{\infty }(4\frac{{\left(-1\right)}^n}{{\left(1+2n\right)}^2}\left(1-\sqrt{{}_1^{2}D}\right)\sin \left[\left(\frac{\pi }{2}+\pi n\right)z\right]\exp \left(-{\left(\frac{\pi }{2}+\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right)+}\\ +\frac{{\left(-1\right)}^n}{n^2}\left(1+\sqrt{{}_1^{2}D}\right)\cos \left(\pi nz\right)\exp \left(-{\left(\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right)).\end{array}$$E179

Returning to the original, we obtain an expression for the first coefficient of expansion of the boundary-layer functions in the layer

$$\begin{array}{l}{\Pi }^{\left(1\right)}=\frac{{}_1^{0}D}{{\pi }^2}\left(\frac{\exp \left(-t\cdot \mathrm{At}\right)}{\sqrt{\pi t}}+\sqrt{\mathrm{At}}erf\sqrt{t\cdot Ad}\right)\times \{4\left(1-\sqrt{{}_1^{2}D}\right)\sin \left(\frac{\pi }{2}z\right)\exp \left(-{\left(\frac{\pi }{2}\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right)+\\ +{\displaystyle \sum _{n=1}^{\infty }(4\frac{{\left(-1\right)}^n}{{\left(1+2n\right)}^2}\left(1-\sqrt{{}_1^{2}D}\right)\sin \left[\left(\frac{\pi }{2}+\pi n\right)z\right]\exp \left(-{\left(\frac{\pi }{2}+\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right)+}\\ +\frac{{\left(-1\right)}^n}{n^2}\left(1+\sqrt{{}_1^{2}D}\right)\cos \left(\pi nz\right)\exp \left(-{\left(\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right))\}.\end{array}$$E180

Equation (180) allows determining the boundary values of the first boundary layer coefficients for the surrounding half-spaces, the problem to determine which are the following:

$$\frac{{\partial }^2{\Pi }_1^{\left(1\right)u}}{\partial z^2}-\left(p+\mathrm{At}\right){\Pi }_1^{\left(1\right)u}=0,z>1,{\Pi }_1^{\left(1\right)u}|{}_{z=1}={\Pi }^{\left(1\right)u}|{}_{z=1},{\Pi }_1^{\left(1\right)u}|{}_{y+z\to \infty }=0,$$E181

$$\frac{{\partial }^2{\Pi }_2^{\left(1\right)u}}{\partial z^2}-{}_2^{1}D\left(p+\mathrm{At}\right){\Pi }_2^{\left(1\right)u}=0,z<-1,{\Pi }_2^{\left(1\right)u}|{}_{z=-1}={\Pi }^{\left(1\right)u}|{}_{z=-1},{{\Pi }_2^{\left(1\right)u}|}_{y+\left|z\right|\to +\infty }=0.$$E182

Solutions to (181), (182) is represented through $${\Pi }^{\left(1\right)и}$$ by the following form:

$${\Pi }_1^{\left(1\right)u}={{\Pi }_{}^{\left(1\right)u}|}_{\text{z}=1}\text{exp}\left(-\sqrt{p+\text{At}}\left(\text{z}-1\right)\right)$$E183

$${\Pi }_2^{\left(1\right)u}={{\Pi }_{}^{\left(1\right)u}|}_{\text{z}=-1}\text{exp}\left(\sqrt{{}_2^{1}D\left(p+\text{At}\right)}\left(\text{z}+1\right)\right)$$E184

which in the originals are represented as

$$\begin{array}{l}{\Pi }_1^{\left(1\right)}=\{\frac{1}{\sqrt{\pi t}}\exp \left(-t\cdot \mathrm{At}-\frac{{\left(z-1\right)}^2}{4t}\right)-\frac{\sqrt{\mathrm{At}}}{2}[\exp \left[\left(z-1\right)\sqrt{\mathrm{At}}\right]\mathrm{erfc}\left(\frac{z-1}{2\sqrt{t}}+\sqrt{t\cdot \mathrm{At}}\right)-\\ -\exp \left[-\left(z-1\right)\sqrt{\mathrm{At}}\right]\mathrm{erfc}\left(\frac{z-1}{2\sqrt{t}}-\sqrt{t\cdot \mathrm{At}}\right)]\}(\frac{4{}_1^{0}D}{{\pi }^2}\left(1-\sqrt{{}_1^{2}D}\right)\exp \left[-{\left(\frac{\pi }{2}\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right]+\\ +\frac{{}_1^{0}D}{{\pi }^2}{\displaystyle \sum _{n=1}^{\infty }(\frac{4{\left(-1\right)}^{2n}}{{\left(1+2n\right)}^2}\left(1-\sqrt{{}_1^{2}D}\right)\exp \left[-{\left(\frac{\pi }{2}+\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right]+}\\ +\frac{{\left(-1\right)}^{2n}}{n^2}\left(1+\sqrt{{}_1^{2}D}\right)\exp \left[-{\left(\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right])),\\ \end{array}$$E185

$$\begin{array}{l}{\Pi }_2^{\left(1\right)}=\{\frac{1}{\sqrt{\pi t}}\exp \left(-t\cdot Ad-\frac{{}_1^{2}D{\left(z+1\right)}^2}{4t}\right)-\\ -\frac{\sqrt{\mathrm{At}}}{2}[\exp \left[-\left(z+1\right)\sqrt{{}_1^{2}D\mathrm{At}}\right]\mathrm{erfc}\left(\sqrt{t\cdot \mathrm{At}}-\frac{\left(z+1\right)\sqrt{{}_1^{2}D}}{2\sqrt{t}}\right)-\\ -\exp \left[\left(z+1\right)\sqrt{{}_1^{2}D\mathrm{At}}\right]\mathrm{erfc}\left(-\frac{\left(z+1\right)\sqrt{{}_1^{2}D}}{2\sqrt{t}}-\sqrt{t\cdot \mathrm{At}}\right)]\}\times \\ \times (\frac{4{}_1^{0}D}{{\pi }^2}\left(\sqrt{{}_1^{2}D}-1\right)\exp \left[-{\left(\frac{\pi }{2}\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right]+\\ +\frac{{}_1^{0}D}{{\pi }^2}{\displaystyle \sum _{n=1}^{\infty }(\frac{4{\left(-1\right)}^{2n+1}}{{\left(1+2n\right)}^2}\left(1-\sqrt{{}_1^{2}D}\right)\exp \left[-{\left(\frac{\pi }{2}+\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right]+}\\ +\frac{{\left(-1\right)}^{2n}}{n^2}\left(1+\sqrt{{}_1^{2}D}\right)\exp \left[-{\left(\pi n\right)}^2\frac{{}_1^{0}D}{\mathrm{Pd}}y\right])).\end{array}$$E186

The solution to the nonstationary problem (44) - (50) in the asymptotic approximation, taking into account the boundary-layer functions, represented as the sum of (134), (161) and (180) in the layer and, relevantly, for the surrounding area

$$\rho ={\rho }^{\left(0\right)}+\epsilon \left({\rho }^{\left(1\right)}+{\Pi }^{\left(1\right)}\right),{\rho }_1={\rho }_1^{\left(0\right)}+\epsilon \left({\rho }_1^{\left(1\right)}+{\Pi }_1^{\left(1\right)}\right),{\rho }_2={\rho }_2^{\left(0\right)}+\epsilon \left({\rho }_2^{\left(1\right)}+{\Pi }_2^{\left(1\right)}\right).$$E187

6. Conclusion

Figure 2 shows the dependence (134), the density of radioactive contaminants on the radial coordinate in the zero approximation for different times of injection with (curves 1, 2, 3) and without (curve 4), radioactive decay, as well as in diffusional approximation (curve 5). In many technological and non-technological liquid wastes of atomic industry, as one of the component contains a radioactive isotope of strontium $$S{r}_{90}^{38}$$ with a half-life $$T_{1/2}=28$$years, which is a very dangerous chemical element for the biological structures due to its ability to replace atoms of calcium. Calculated options: half-thickness of the layer h = 1m, the diffusion coefficients D_1z= D_2z =10^-11m²/c, D_z=10^-9 m²/c; borehall radius of 0.1 m, the polluter - Strontium 90; activity of the solution - 1 Ci / l (high level waste ), the initial density of radioactive substances in the solution 7.34 g/m³; injection volume 100 m³/day (At = 75.8, Pe =1.84 10⁷)

Figure 2.

Plots of density on distance for different times of observation: Given the radioactive decay of 1 - 5 years, 2 - 15 years old, 3 - 25 years, 4 - 25 years (in the absence of decay At = 0), 5 - 25 years (in the absence of decay in non-diffusive approximation)

Figure 2 allows carrying out estimate of the contribution of diffusion and radioactive decay. The figure shows that the account of the radioactive decay (see curves 3 and 4) is the priority in comparison with the account of the mass transfer of the layer (curves 4 and 5) with the surrounding formations by diffusion. The figure also implies that the contribution of diffusion is essential in the front zone of radioactive contamination, where it is comparable with the magnitude of the density of matter in solution. With the zero approximation the size of the zone of contamination are determined.

Figure 3.

The dependence of the density of radioactive contaminants on the radial coordinate for the dimensionless time t = 0.01 (a), t = 0.001 (b): 1.2 - without and with taking into account the boundary-layer correction, respectively, 3 - zero approximation at z = 0

On the fig. 3 a, b it is shown that taking into account the boundary layer solution eliminates the disadvantage of the first approximation, which consists in the fact that the density of radioactive contaminants exceeds one (curve 1). In the calculations agreed: $$\text{Pd}=2\cdot {10}^7$$, $${}_0^{1}D=0.01$$, $${}_1^{2}D=1$$, $$\epsilon =0.05$$, Ad = 2200, which corresponds to the half-life $$T_{1/2}=1\text{year}$$. Boundary layer correction results in the behavior of the curve in accordance with the conditions of the problem (see curves 1, 2) and at the same time clarifies the first approximation. Note that for short times accounting of the first expansion coefficient leads to refinements of the design parameters up to 20% (curves 1 and 3 in Fig.3 b). At the same time the importance of taking into account the boundary-layer functions in the near-well zone is illustrated ($$r<\text{70}$$).

As the figure 3 b shows, on the border of the front distribution of contaminants occurs a sharp jump of density, indicating the presence of internal moving boundary layer, the patterns of which have not been investigated, and mathematical methods of eliminating it - not developed.

Appendix: Note on the stationary solution to the problem

It is extremely important to find a stationary solution to set the maximum size of the zone of contamination. The equations describing the stationary regime are obtained from (44) - (50) if we consider the first term $$\left(\partial {\rho }_i/\partial t\right)$$ to be zero. The solution to this problem is given in [39, 40, 43]. Here we note that the solution to the stationary problem can be obtained from the first asymptotic approximation for $$t\to \infty$$.

Figure 4.

The dependence of the densities of radioactive contaminants in the porous layer for the stationary case (zero approximation) on the distance from the borehole at different decay constants: 1 - Ad = 0.01, 2 - 0.1, 3 - 1. Other parameters: Pd = 10², $$\delta ={10}^{-4}$$, $$D_1^2=1$$

Fig. 4 shows the dependence of the density of radioactive contaminants on the radial coordinate at the center of layer for various radioactive contaminants: curve 1 - ²³⁹Pu (T_1/2 = 2,24 10⁴ years), curve 2 – ²²⁶Ra (T_1/2 = 1590 years), curve 3 – ⁹⁰Sr (T_1/2 = 28 years).

The zero approximation in this case is the most important; it determines the general form of the dependence. The value of the density of the pollutant decreases exponentially, and as follows from the graphs, even for the middle half-life and most dangerous radionuclides (⁹⁰Sr, ¹³⁷Cs) at distances of 200 h (200 m) of the order of percent of maximum, observed in the area of injection.

Figure 5.

The dependence of the density radioactive contaminants in the stationary case (zero approximation) on the z coordinate at different distances from the borehole: 1 - r = 0, 2 – 100, 3 - 200. Other parameters: Ad = 0.01, Pd = 10², $$\delta ={10}^{-4}$$, $$D_1^2=1$$

Figure 5 shows a picture of the field distribution of the radioactive contaminant in the stationary case along the vertical coordinate (zero approximation). "Slices" are given for distances 0, 100h and 200h from the axis of borehall. We see that for the middle half-life nuclides (Т_1/2 ~ 30 years) in the covering and underlying layers pollutant density decreases rapidly, and even at distances of 0,5h are negligible.

In general, the increase in the parameter Pd (rate of the injection of the solutions) leads to the "elongation" of the graphic along the radial direction, reducing the At (which corresponds to an increase in the average lifetime of the nuclide) - to "enhance" the graph along the axes r and z. The field of pollutant remains limited.

Figure 6.

The dependence of the density of radioactive contaminant on the distance to the borehole axis Graphs are constructed (for the dimensionless time t = 100): the grid method - 1 and the method of asymptotic expansion - 2. Other parameters: At = 0.1, Pd = 102, $$\delta ={10}^{-3}$$, $$D_1^2=1$$

Fig. 6 shows the results obtained using a modified method of asymptotic expansions and the numerical solution to the problem of mass transfer by the grid method. Numerically was solved the problem (67) - (74), neglecting the radial diffusion.

Comparisons of the curves shown in Fig. 6 allow to conclude that the results obtained by numerical and asymptotic methods are in a good agreement.

So, based on the asymptotic method, approximated analytical solution to the problem of subterranean waste disposal is obtained, and accounting the boundary layer correction allows to provide the calculation of the areas of radioactive contaminants in the subterranean horizons with high accuracy at all distances from the injection borehall, and thus to clarify the forecast of the areas of radioactive contamination to ensure the environmental safety.

In conclusion, note that the above modification of the asymptotic method is quite general and provides the construction of "exact on the average" analytical solutions as to the nonstationary problem of the underground waste disposal as well as to the other problems of underground thermo- and hydrodynamics. The zero approximation of the asymptotic solution is of the particular importance, because it describes the average value of the variables, which is important for many practical problems.

List of designations

A, B, C, E, F,M	–	auxiliary functions;
$$\alpha$$, Ad, At	–	dimensional and dimensionless constants for radioactive decay for the diffusion and the temperature problem, respectively, 1/c;
$${\lambda }_{z1},{\lambda }_z,{\lambda }_{z2}$$,; $${\lambda }_{r1},{\lambda }_r,{\lambda }_{r2}$$		coefficients of thermal conduction for the covering, porous, and the underlying layers in the vertical and radial directions, respectively, W/(m·K);
$${\rho }_{\text{п}}$$, $${\rho }_{\text{п1}}$$, $${\rho }_{\text{п2}}$$	–	density of the porous, covering and underlying layers, kg/m^3;
$${\rho }_{\text{d}}$$, $${\rho }_{1\text{d}}$$, $${\rho }_{2\text{d}}$$	–	dimensional concentration of the impurity in a porous, covering and underlying layers, kg/m^3;
$${\rho }_f$$, $${\rho }_s$$		dimensionless densities of the impurity in the carrier, the skeleton;
Dr,D1r, D2 ,; Dz,D1z, D2z	–	coefficients of diffusion of the layers in the radial and vertical directions, m^2/s;
δi,j	–	the Kronecker delta,
$${\text{μ}}_s,{\text{μ}}_w$$	–	chemical potentials of the skeleton and water, respectively
$$g\left({\text{μ}}_s,{\text{μ}}_w\right)$$	–	function of the mass transfer between the skeleton and the fluid;
H, h	–	power and half-thickness of the porous layer, m;
KГ	–	Henry’s coefficient;
$$\stackrel{⌢}{L}$$	–	differential operator;
Pd	–	analog of the parameter Peclet;
$$R_p$$	–	radiuses of the radionuclide contamination;
p, s	–	parameters of the Laplace-Carson,;
r0	–	radius of the bonehole, m;
rd, zd, r, z	–	dimensional and dimensionless cylindrical coordinates, m;
$$\tau$$, t	–	dimensional and dimensionless time, s;
m, m0 , m1	–	effective, initial and maximum porosity;
$$\text{v}$$	–	filtration rate, m/s;
$$\text{<msup>v <mo>′</mo></msup>}$$	–	rate of convective transport of contaminant in the porous layer, m/s;
w	–	the true velocity of the fluid, m/s;
qd, q	–	dimensional and dimensionless source function (mass, kg/(s·m^3);)
$$k,k_1,k_2$$	–	stoichiometric coefficients;
$$\theta \text{,}{\theta }_{\text{1}},{\theta }_{\text{2}}$$	–	the remainder term of the asymptotic expansion in a porous, covering and underlying layers.