Estimation of Radon Flux Density Changes Using Mathematical α-Model Before and After a Strong Earthquake in Kamchatka

1. Introduction

In today’s world, one of the most dangerous threats to human life is seismic activity, the effects of which can be extremely destructive.

In the Russian segment of the Western Pacific continent-ocean transition zone, two seismogenic structures can be identified that differ significantly in energy and spatial manifestation of seismicity, the Kurilo-Kamchatka and Sakhalin seismic zones. While the former is a subduction zone with possible strong earthquakes with magnitudes up to $M_w=9.0$ occurring at depths of 0-700 km, the latter is confined to the western boundary of the Okhotsk lithospheric microplate and is characterized by crustal earthquakes.

For the Kamchatka region the problem of seismicity has the highest priority, as the Kamchatka peninsula is located in one of the most earthquake-prone areas of the Earth. The cities of Petropavlovsk-Kamchatsky, Yelizovo and Vilyuchinsk are located in the zone of 9-point earthquakes, which once again proves the need to develop effective methods of seismic activity forecasting. Strong seismic events can cause a tsunami wave, which can cause great economic damage to the shores of Kamchatka and the Kuril Islands.

A possible way to solve this problem, the development of an effective forecasting system based on an in-depth understanding of the processes occurring in the Earth’s crust in preparation for a future earthquake focal.

One of the methods of studying such processes is monitoring changes in the concentration of radon gas (${}^{\mathtt{222}}\mathtt{Rn}$) in the subsurface, atmospheric air and water [1]. Such changes can serve as precursor signals of preparation of strong earthquakes. The study of ${}^{\mathtt{222}}\mathtt{Rn}$ variations in subsurface air helps to obtain information about geodynamics processes in the mountain massive and properties of the medium through which subsurface gases migrate. Long-term observations and their comparison with the seismicity of the region allow tracking changes in the character of geodynamics processes. This information can be used to search for precursor anomalies that may indicate the approach of a strong earthquake.

One of the uranium decay products is the noble gas ${}^{\mathtt{222}}\mathtt{Rn}$. The decay of radon and its products occurs with the emission of three types of radiation, by which its concentration can be measured. Radon is transported into the atmosphere by diffusion and convection.

The Petropavlovsk-Kamchatka geodynamics test site is monitoring ${}^{\mathtt{222}}\mathtt{Rn}$ on a network of points. The network is created in such a way that the main attention is directed to the Avacha and Kronotsky Bays and the south of Kamchatka. The observation points are equipped with complexes for recording subsurface gases, primarily radon and hydrogen. The network has been operating in different configurations for more than 20 years. On the basis of the obtained data, including volumetric radon activity (RVA) for a number of earthquakes in the Avacha Gulf region, two types of precursor anomalies were identified in the subsurface radon field.

Type A anomalies are characterized by their in-phase manifestation at 3-5 points with a relative shift in time. The most common forms of these anomalies are bay-shaped and stepped, of different polarity. The precursor anomaly of type A was also registered for the deep (origin depth 177 km) Zhupanovskoe earthquake with $M_w=7.2$ and for the earthquake on 03 April 2023 with $M_w=6.9$. The possible reason for the appearance of in-phase anomalies on the network of points of registration of subsurface gases are deformation processes such as a solitary ‘deformation wave’, which can arise due to quasi-viscous flow of geomaterial at the last stage of earthquake preparation.

Type B precursor anomalies were registered, as a rule, in one point of the network, had a definite shape and stood out well against the general background. Calculations performed by mathematical modeling methods showed that the mechanism of occurrence of such anomalies is associated with the injection of ${}^{\mathtt{222}}\mathtt{Rn}$ into the groundwater flow with full transverse mixing under the influence of a deformation stress pulse [2].

Identification of signs of such processes in the data of radon observations is one of the actual tasks of geophysics. Application of new experimental data of subsurface gas variations in comparison with data on seismicity of the region and study of ${}^{\mathtt{222}}\mathtt{Rn}$ migration processes with the help of the most modern methods of mathematical modeling in order to interpretation of anomalous effects preceding earthquakes is a new and actual method of research in Earth sciences and in particular in development of earthquake forecasting methods.

The earthquake with magnitude $M_w=7.0$ has happened on August 17, 2024 at 19:10:26 (UTC) in the northern part of Avacha Gulf, not far from Shipunsky peninsula, later named Shipunsky in connection with geographical location of epicentre (52.931°N, 160.133°E). The epicentral distance to Petropavlovsk-Kamchatsky was 100 km and the depth was 29 km according to the National Earthquake Information Centre (NEIC, U.S. Geological Survey). The intensity of shaking in Petropavlovsk-Kamchatsky was VI-VII points according to seismic intensity scale (SIS).

Before this earthquake, a number of anomalous effects were registered in the subsurface radon field at the network of points [3] near Petropavlovsk-Kamchatsky. The data of radon monitoring before and after this earthquake were used in this work to estimate the changes in the radon flux density in the accumulation chamber.

On (Figure 1) shows the RVA and atmospheric pressure curves. On the RVA curve for estimation of radon flux density, two least noisy areas were selected for exogenous impacts. In the second plot (Figure 1 b), the abrupt release of RVA is associated with the arrival of a cyclone and a decrease in atmospheric pressure. To use these data for modeling purposes, the RVA data were removed at the specified site and the skip was filled in using the Lubushin A.A. software package. [4].

The figure shows that three days before the earthquake there was a sharp decrease in RVA, and then followed by a new stage of accumulation. It is supposed that such behavior of RVA is connected with the last stage of development of the source of future earthquake and is determined by propagation of deformations in the Earth crust. Similar behavior of RVA curves was recorded before a number of strong earthquakes in southern Kamchatka and earlier [3].

The article has the following structure: section 1 presents the introduction, describes the scope and object of the study, and presents the structure of the article; section 2 briefly presents information about radon monitoring as the main method of obtaining experimental data on RVA accumulation used in the study; section 3 describes concepts and models used in mathematical modeling of the RVA accumulation process; section 4 presents the classical mathematical model of RVA, as well as concepts and quantities related to this process, describes the hereditary $\alpha$-model RVA as a generalization of the classical model, formulates the direct problem and the method of its solution; section 5 formulates the inverse problem of restoring the values of several parameters ${\lambda }_0$ and $\alpha$ and describes the method of its solution; section 6 describes the method of estimating the flux density ${}^{\mathtt{222}}\mathtt{Rn}$ based on the used mathematical models and their parameters; section 7 presents the results of solving inverse problems for ${\lambda }_0$ and $\alpha$ on the basis of different experimental data of RVA, and estimates the flux density ${}^{\mathtt{222}}\mathtt{Rn}$ on the basis of the recovered values of the hereditary $\alpha$-model RVA; section 8 describes the research contribution and the various sources of funding; section 9 summarizes the study, formulates the conclusions.

2. About Radon Monitoring

It is known that radioactive metal ${}^{\mathtt{226}}\mathtt{Ra}$ is constantly present in the Earth’s crust. The half-life of this element produces ${}^{\mathtt{222}}\mathtt{Rn}$ (radon), an inert radioactive gas with a half-life of about 3.85 days. This led to the development of radon monitoring, a special area of emanation and radiation monitoring of geo-medium and processes. Radon monitoring includes measurement of the level of ${}^{\mathtt{222}}\mathtt{Rn}$ and other radon-forming elements, as well as their distribution in space.

On the Kamchatka Peninsula, a network of observation sites [5] has been deployed to monitor ${}^{\mathtt{222}}\mathtt{Rn}$ and additional environmental parameters (temperature, pressure, moisture, etc.) that can affect the resulting accumulation curve using accumulation chambers with gas-discharge counters. When planning the locations of subsurface gas monitoring network stations in Kamchatka, special attention is paid to river valleys that run along crustal faults. These zones are characterized by increased permeability, which creates favorable conditions for subsurface gases to escape into the atmosphere, as described in the article [6].

Remark 1.

Analysis of RVA data and related parameters obtained during continuous monitoring is one of the methods of searching for earthquake precursors.

This is due to the fact that RVA is affected by changes in the stress-strain state of the medium through which subsurface gas comes to the surface [7], so ${}^{\mathtt{222}}\mathtt{Rn}$ is considered a well-known and well-proven indicator of the processes occurring in such a medium [8,9]. Monitoring of ${}^{\mathtt{222}}\mathtt{Rn}$ as a method of searching for precursors of seismic events has proved itself in recent years [10,11], especially as a short-term precursor (up to 15 days) [12,13].

All points in the network are organized according to a common principle. The sensors are usually located in accumulation chambers (10 liter galvanized bucket) in the aeration zone. When the equilibrium between ${}^{\mathtt{222}}\mathtt{Rn}$ and its decay products is reached in the accumulation chamber, the intensity of $\beta$-radiation increases, which increases the sensitivity of measurements. The insignificant convective component of subsurface air in the chamber was ensured by holes in the bottom of the bucket and loose packing of the pipe through which the detector was lowered into the chamber (Figure 2)

Gas-discharge counters are the most common detectors of $\beta$- and $\gamma$- radiation. The high sensitivity of counters allows to register single quanta of ionizing radiation, and the large output signal is easily registered by recalculation circuits. All this allows passive registration of ${}^{\mathtt{222}}\mathtt{Rn}$ in subsurface air by $\beta$-radiation of short-lived products of its decay with a high degree of reliability and rather simple metrology.

Conversion from pulses to RVA is carried out by the empirical formula $A\left(t\right)=9·M$ (Bq/m${}^3$), where M is the number of pulses registered by $\beta$-radiation sensors per minute), obtained as a result of synchronous registration by certified radiometers RS-410F by femto-TECH (USA), RRA-01M-03 (Russia) and recorder with GDC-19 detector.

The data used in this work were obtained from the RVA sensor at the INSR point installed in the aeration zone of loose sediments at a depth of 3 m from the day surface [3].

Figure 3. Experimental RVA data for 14 days, from the INSR observation point obtained before the earthquake during the period: 27 July 2024 (16:00) – 10 August 2024 (15:30).

Figure 3. Experimental RVA data for 14 days, from the INSR observation point obtained before the earthquake during the period: 27 July 2024 (16:00) – 10 August 2024 (15:30).

Figure 4. (green) – experimental RVA data for 13 days, from INSR observation point obtained after the earthquake (include earthquake) during the period: 16 August 2024 (16:00) – 29 August 2024 (15:30); (red) – 7.0 magnitude earthquake.

Figure 4. (green) – experimental RVA data for 13 days, from INSR observation point obtained after the earthquake (include earthquake) during the period: 16 August 2024 (16:00) – 29 August 2024 (15:30); (red) – 7.0 magnitude earthquake.

3. Approaches used for RVA modeling

The mechanism of ${}^{\mathtt{222}}\mathtt{Rn}$ transport in geological media is usually considered within the framework of the [14] emanation method and many physical and mathematical models of the [15,16,17,18].

We consider the classical mathematical model of RVA in time $A\left(t\right)$ in an accumulation chamber based on the equation with ordinary derivative of the 1st order (ODE-model RVA), as well as concepts and quantities related to this process [3,19,20]. This model is the basis for a more generalized process model developed by the authors[2,21]. The weakness of the ODE-model RVA is that it does not take into account changes in the stress-strain state of the medium, assuming that the transfer process ${}^{\mathtt{222}}\mathtt{Rn}$ occurs in a homogeneous medium.

Remark 2.

In the simplest case, the RVA ODE model assumes a stationary accumulation mode in the chamber, which means that the model must take into account at least ${\lambda }_0$ – a constant coefficient of air exchange ratio (AER). This is one of the main parameters of both considered RVA models, both ODE and its generalization.

Actually, in order to take into account that the transport of ${}^{\mathtt{222}}\mathtt{Rn}$ takes place in an heterogeneous medium, the authors earlier [2,21] proposed generalizations of the ODE model RVA using the fractional derivative (FD) apparatus [22,23]. The generalization to the hereditary $\alpha$-model RVA consists in replacing the 1st order ordinary derivative by the Caputo (Gerasimov-Caputo) FD [24,25] of constant $0<\alpha <1$ order.

The prerequisites for this are several points. First, under the action of tectonic tensions there is a change in the vertical velocity of gas flow[26]. Secondly, the stress-strain state of the medium leads to changes in such characteristics of the geological medium as porosity, permeability, fracturing [2], and the transfer process ${}^{\mathtt{222}}\mathtt{Rn}$ occurs in such a permeable medium [27].

Remark 3.

Therefore, the study uses an hereditary α-model RVA to account for the variation of the transport intensity ${}^{\mathtt{222}}\mathtt{Rn}$ during the accumulation process developed during the [2,21,28] research cycle.

Using data when the RVA accumulation curve is clearly visualize and free of noise will allow us to construct model curves that are as close as possible in shape to the real data and with maximum correlation. Therefore, for modeling from (Figure 1) data were selected from time intervals when there were no abrupt changes in atmospheric pressure or winds from cyclones passing by the observation point. Both factors cause blowing of ${}^{\mathtt{222}}\mathtt{Rn}$ out of the chamber, due to which the RVA curves show <<bumpy>> regions that will interfere with the work of mathematical modeling algorithms, reduce the degree of correlation between model curves and data, and also reduce the accuracy of reconstructed values of $\alpha$ and ${\lambda }_0$.

The restoration of the parameters $\alpha$ and ${\lambda }_0$ is carried out by methods for solving inverse problems [29,30] from known experimental RVA accumulation data before (Figure 3) and after (Figure 4) earthquake in order to recover the values of the constant parameters ${\lambda }_0$ and $\alpha$ in the hereditary $\alpha$-model RVA, thus testing the hypothesis about the influence of the change in $\alpha$ – transport intensity due to the seismic event.

4. RVA models and direct problem statement

In our study, we depend on a well-studied mathematical model of ${}^{\mathtt{222}}\mathtt{Rn}$ accumulation in an accumulation chamber that effectively describes the accumulation modes [3,19,20]. This model is a linear ordinary differential equation (ODE) of the following form:

$${\displaystyle \frac{dA\left(t\right)}{dt}}=S\left(t\right)-\left(A\left(t\right)-A_{Rn}^{atm}\right){\lambda }_v\left(t\right)-\left(A\left(t\right)-A_{Rn}^{atm}\right){\lambda }_{Rn},$$

(1)

where

• $A\left(t\right)\in C^1[0,T]$ – time dependence of RVA in the chamber, [Bq/m${}^3$];
• $C^1[0,T]$ – class of continuous-differentiable functions;
• $A_{Rn}^{atm}$ – RVA in outside air, [Bq/m${}^3$];
• $t\in [0,T]$ – the current RVA simulation time, [s];
• $T>0$ – total simulation time, [s];
• $S\left(t\right)$ – some function that describes the rate of entry ${}^{\mathtt{222}}\mathtt{Rn}$, i.e., the total specific entry per unit volume of the chamber, [Bq/m${}^3$s];
• ${\lambda }_v\left(t\right)$ – a function describing the time dependence AER in the chamber, [s${}^{-1}$];
• ${\lambda }_{Rn}=2.1·{10}^{-6}$ – radon decay constant, [s${}^{-1}$].

As shown in [3,20], the (1) model can be significantly simplified in certain situations. For example, the term at ${\lambda }_{Rn}$ of the equation (1) can be ignored, since even a fully enclosed room with sensors has ${\lambda }_v\left(t\right)\ge 0.05$ [h] (i.e., $\approx 1.3889·{10}^{-5}$ [s]), which is at least an order of magnitude larger than ${\lambda }_{Rn}$, and thus the effect of the 3rd term of (1) on the calculations of this model is small.

Remark 4.

It is known that the transport of ${}^{\mathtt{222}}\mathtt{Rn}$ in the vertical direction may be due to thermal-fluid convection, turbulent effects due to changes in meteorological factors, diffusion due to the pressure gradient in the Earth’s crust, diffusion due to the concentration gradient of ${}^{\mathtt{222}}\mathtt{Rn}$, and others. [31]. Therefore, in the (1) model, the $S\left(t\right)$ quantity describes simultaneously the convective $S^C\left(t\right)$ and diffusive $S^{\Delta }\left(t\right)$ mechanisms. Since there is no significant convective flow of subsurface air from the surface below the accumulation chamber, the term $S^C\left(t\right)$ can be neglected. Which mainly characterizes the accumulation regimes in accumulation chambers.

Remark 5.

The modeled accumulation process ${}^{\mathtt{222}}\mathtt{Rn}$ can be considered stationary in the sense of the arrival rate ${}^{\mathtt{222}}\mathtt{Rn}$ when the radon flux density (RFD) from the surface under the accumulation chamber is constant and when there are no dramatic changes in AER. Then the model parameters ${\lambda }_v\left(t\right)={\lambda }_0$ and $S\left(t\right)=S$ are constant values, and the RVA will have an accumulative character with saturation yield: $A_{max}=S/{\lambda }_0$ [Bq/m${}^3$] [3]. From where we obtain that $S=A_{max}{\lambda }_0$.

Considering the above, the Cauchy problem for (1) takes the form:

$${\displaystyle \frac{dA\left(t\right)}{dt}}=-{\lambda }_{0}A\left(t\right)+A_{max}{\lambda }_0,A\left(0\right)=A_0,$$

where $A_0$ – a known constant, the value of RVA at time $t=0$.

Remark 6.

Further, when working with experimental data, it will be necessary to make for the time series RVA normalization to the maximum, thus passing for the notation $A\left(t\right)$ to <<relative units>> [Rel.unit].

Then $A_{max}=1$, and hence the ODE-model will take the form:

$${\displaystyle \frac{dA\left(t\right)}{dt}}={\lambda }_0\left(1-A\left(t\right)\right),A\left(0\right)=A_0,$$

(2)

Remark 7.

Further, the stationary classical RVA model (2) will be shortly named the ODE-model RVA.

For the Cauchy problem (2), an analytical solution can be obtained:

$$A\left(t\right)=1-e^{-{\lambda }_{0}t}\left(1-A_0\right),A\left(0\right)=A_0.$$

(3)

where at $t\to \infty$ we can see that the RVA $A\left(t\right)\to 1$.

The (1) and (2) equations use ODE, which significantly limits the flexibility of the ODE model. Therefore, the authors in [2,21] propose a modification of the (2) model consisting in replacing the 1st order ordinary derivative by a fractional derivative of [22,23] of constant order.

Remark 8.

It is assumed that porosity of the medium is due to the presence of isolated pores, and permeability of the medium is understood as the presence of channels conducting gas between the pores. The porosity of the medium can lead to slowing down of the gas transfer process, i.e., subdiffusion, while the permeability of the medium, on the contrary, leads to acceleration, i.e., superdiffusion [14]. Such processes belong to the phenomena of anomalous diffusion [32].

Anomalous diffusion can be related to the property of a system or a medium to remember for some time the impact made on it - [33]. Therefore, the fractional derivative Caputo (Gerasimov-Caputo) [24,25] of constant order is introduced in the mathematical model RVA to describe the heredity:

$${\partial }_{0,t}^{\alpha }A\left(\phi \right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{\displaystyle \frac{dA\left(\phi \right)}{dt}}{\displaystyle \frac{1}{{(t-\phi )}^{\alpha }}}d\phi ,0<\alpha <1,$$

(4)

where $\Gamma (·)$ is the Euler gamma function:

$$\Gamma \left(x\right)={\int }_0^{\infty }{e}^{-t}{y}^{x-1}dy,x\in \mathbb{C}:Re\left(x\right)>0.$$

Remark 9.

There are other definitions of the fractional order derivative. They can be found, for example, in scientific works [34,35]. However, the use of fractional Gerasimov-Caputo derivative (4) is the simplest and at the same time quite effective approach.

A series of works by the authors of [2,21,28] is devoted to investigating issues related to mathematical modeling of RVA, where it is assumed that the parameter $\alpha$ describes the fractality of the [36] of the geo-medium and is related to its characteristics such as porosity, permeability and fracturing. The hereditary $\alpha$-model RVA is represented by the following Cauchy problem:

$${\varphi }^{1-\alpha }{\partial }_{0,t}^{\alpha }A\left(\phi \right)=-{\lambda }_{0}A\left(t\right)+{\lambda }_0,A\left(0\right)=A_0,$$

(5)

where the difference from (1) is important to note:

• $A\left(t\right)\in C^2[0,T]$ – the solve function, the dependence of RVA on time in the chamber, [Rel.unit];
• $\alpha$ – transport intensity ${}^{\mathtt{222}}\mathtt{Rn}$, order of fractional derivative (4), dimensionless constant;
• $C^2[0,T]=\mathbb{A}$ – a class of twice continuously differentiable functions;
• $\varphi$ is some positive constant with time dimension [37]. In the following, the work is carried out with undimensionazed and normalized to the maximum, so we assume $\varphi =1$.

Definition 1.

It is assumed that the order of the fractional derivative α - responsible for the intensity of the transport process ${}^{\mathtt{222}}\mathtt{Rn}$, is related to the characteristics of the geosphere: permeability, porosity, fracturing, and possibly also to the fractal dimension, etc. It is assumed that the constant $0<\alpha <1$ describes the average permeability value of the geosphere under the accumulation chamber of the recording RVA over some volume. Moreover, if $\alpha \to 0$, it characterizes less permeable for ${}^{\mathtt{222}}\mathtt{Rn}$ medium, and if $\alpha \to 1$, it characterizes more permeable medium.

Remark 10.

The hereditary α-model RVA (5) solved by the scheme (7) at the value $\alpha =1$ will pass to the ODE-model RVA (3), this is demonstrated in [28,38], which suggests that the generalization of (5) is correct.

To numerically solve the problem (5), we will use the previously developed nonlocal implicit finite-difference scheme (IFDS) [39] defined in a uniform grid domain:

$$\begin{array}{c}\hfill \begin{array}{c}h=T/N,\widehat{\Omega }=\left\{(t_i=ih):0\le i<N\right\},\widehat{\mathbb{A}}\in \widehat{\Omega },\\ A\left(t\right)=A_i,0<A_i<1.\end{array}\end{array}$$

(6)

Definition 2.

Then the difference direct problem:

$$\begin{array}{c}\hfill \begin{array}{c}{A}_i=1-{\displaystyle \frac{\widehat{{\partial }_{0,ih}^{\alpha }}{A}_i}{{\lambda }_0}},A_0-const,1\le i<N,\\ \widehat{{\partial }_{0,ih}^{\alpha }}{A}_i={\displaystyle \frac{h^{-\alpha }}{\Gamma (2-\alpha )}}\sum _{j=0}^{i-1}\left({(j+1)}^{1-\alpha }-j^{1-\alpha }\right)\left(A_{i-j}-A_{i-j-1}\right).\end{array}\end{array}$$

(7)

is a Cauchy problem consisting in finding a discrete function $A_i$ belonging to a known class $\widehat{\mathbb{A}}\in \widehat{\Omega }$, with known constants α and ${\lambda }_0$.

IFDS scheme (7) and its generalized analogs have been tested in a number of test and applied problems [38,39]. Then (7) is solved by the modified Newton method (MNM), and then the IFDS-based difference direct problem is unconditionally stable.

5. Inverse problem on parameters ${\lambda }_0$ and $\alpha$ for the hereditary $\alpha$-model RVA

Previously, in [2,21,28] the parameters of models ${\lambda }_0$ and $\alpha$ were unknown, and were selected by brute force based on some considerations about the process flow, which is a time-consuming approach. The accuracy of the selection was evaluated by the maximum of $R^2$ – the coefficient of determination [40] and $\sigma$ – the Pearson correlation coefficient [41] with the maximum normalized experimental RVA data. This naturally leads to the search for solutions to automate the selection of optimal parameters.

Let $A_i\in \widehat{\mathbb{A}}$ (and respectively $A\left(t\right)\in \mathbb{A}$) is a function of some known class, but its solution depends on the set of parameters $\overrightarrow{X}=\left[X_0,....,X_{K-1}\right]$, where $K=2$, and $X_0=\alpha$, $X_1={\lambda }_0$. Let the values of the discrete decision function $A_i\in \widehat{\mathbb{A}}$ are unknown, but additional information (experimental data of RVA) $A_i={\theta }_i=\overrightarrow{\theta }$ on the solution of the difference direct Cauchy problem (7) for the hereditary $\alpha$-model RVA is known.

Definition 3.

Then the difference inverse problem for (7) is to recover the values $\overrightarrow{X}=\left[X_0,X_1\right]$ from the known $A_i={\theta }_i$ experimental RVA data:

$$\begin{array}{c}{A}_i=1-{\displaystyle \frac{\widehat{{\partial }_{0,ih}^{X_0}}{A}_i}{X_1}},A_i={\theta }_i,1\le i<N,\\ \widehat{{\partial }_{0,ih}^{X_0}}{A}_i={\displaystyle \frac{h^{-X_0}}{\Gamma (2-X_0)}}\sum _{j=0}^{i-1}\left({(j+1)}^{1-X_0}-j^{1-X_0}\right)\left(A_{i-j}-A_{i-j-1}\right).\end{array}$$

(8)

To solve (8), we turn to the unconditional optimization theory [42]. For this purpose, it is necessary to minimize the bias functional:

$$\begin{array}{c}\overrightarrow{\eta }=\overrightarrow{\theta }-\omega \left(\overrightarrow{X}\right),min\left(\Psi \left(\overrightarrow{X}\right)\right)={\displaystyle \frac{1}{2}}\sum _{i=0}^{N-1}{\eta }_i^2={\displaystyle \frac{1}{2}}\sum _{i=0}^{N-1}{\left({\theta }_i-{\omega }_i\right)}^2,\end{array}$$

(9)

where $\overrightarrow{\eta }$ is a bias vector of dimension $N>K$, and the vector $\omega \left(\overrightarrow{X}\right)=\left[{\omega }_0,...,{\omega }_N\right]$ is a vector of model data, i.e. the solution of the difference direct problem (7) with respect to some approximation $\overrightarrow{X}$ obtained by solving the inverse problem.

The difference inverse problem is solved by a Newton-type unconditional optimization method [43], namely the Levenberg-Marquardt iterative method [44,45], represented as:

$$\Delta X=\left(-H^{-1}\right)\times \left(J^T\times \overrightarrow{\eta }\right),H=J^T\times J+\gamma E,$$

(10)

where

• $\Delta X$ – the optimal increment of $\overrightarrow{X}$ for the next iteration;
• E is a unit matrix of dimension $K\times K$;
• $J=J\left(\overrightarrow{X}\right)$ – a Jacobi matrix of dimension $N\times K$ with elements calculated by the formula: $J_{i,k}={\displaystyle \frac{\partial {\eta }_i}{\partial X_k}},i=0..N-1,k=0..K-1$;
• the derivative ${\displaystyle \frac{\partial {\eta }_i}{\partial X_k}}$ is approximated by the difference operator $J_{i,k}={\displaystyle \frac{{\eta }_i^{\delta }-{\eta }_i^{}}{\delta X_k}}$, where $\delta X$ – a given small increment of $\overrightarrow{X}$;
• $\gamma$ is the regularization parameter of the method. If $\gamma \in {\mathbb{R}}_{>0}$ and the Hesse matrix H is positive definite, then $\Delta X$ is the direction of descent for the optimal step of the method;
• Start value: ${\displaystyle {\gamma }^{\left(0\right)}=v·\underset{i}{max}\left(diag\left(J{\left(X^{\left(0\right)}\right)}^T\times J\left(X^{\left(0\right)}\right)\right)\right)}$, where v is a given starting constant.

Remark 11.

The solution of the inverse problem (8) by the Levenberg-Marquardt method (10), hereafter (IP-LB), is reduced to starting from the given constants $X^{\left(0\right)}$, $\delta X$, v, and c – constants for recalculating γ during the loop, by recomputing the solution of the difference direct problem (7) many times with approximations $\overrightarrow{X}$ obtained during the solution of the inverse problem, to compute the optimal values of $\overrightarrow{X}$. More details on the algorithm for implementing the optimization of the vector $\overrightarrow{X}$ can be found in the article [46].

The criterion for obtaining the optimal value is $\epsilon \le \Sigma$, where $\Sigma$ is the given accuracy of IP-LB solution, $\epsilon ={\displaystyle \frac{1}{N}}{\sum }_{i=0}^{N-1}{\left[{\eta }_i^{\Delta }\right]}^2$ – mean square error (MSE) between experimental and model RVA data.

6. About the Flux Density Estimation of ${}^{\mathtt{222}}\mathtt{Rn}$

Variations of RVA are related to such indicator as radon flux density (RFD) from the surface under the accumulation chamber. Taking into account Remarks 4, 5 and according to [3] the arrival rate ${}^{\mathtt{222}}\mathtt{Rn}$ is described by the relation: $S=q\Pi /V$. Then, to estimate RFD based on a model stationary in the sense of the arrival rate ${}^{\mathtt{222}}\mathtt{Rn}$, it suffices to express:

$$q={\displaystyle \frac{A_{max}{\lambda }_{0}V}{\Pi }},$$

(11)

where

• $S=A_{max}{\lambda }_0$ – the arrival rate of ${}^{\mathtt{222}}\mathtt{Rn}$ into the chamber, [Bq/m${}^3$s];
• q – RFD from the surface under the accumulation chamber, [Bq/m${}^2$s];
• $\Pi$ – area of flow under the chamber (Figure 2), [m${}^2$];
• V – volume of the accumulation chamber (Figure 2), [m${}^3$].

Remark 12.

The parameters of the accumulation chamber at the INSR observation point (Figure 2) take the values: $V=0.01$, $\Pi =0.05$.

It is also important to note that when working with experimental data RVA is normalized to the maximum $A_{max}=1$ [Rel.unit]. But if we want to estimate RFD for the original RVA data, we should substitute the maximum un-normalized value in [Bq/$m^3$] in place of $A_{max}$.

However, the values of AER ${\lambda }_0^{\left(1\right)}$ before the earthquake (Figure 3) and AER ${\lambda }_0^{\left(2\right)}$ after the earthquake (Figure 4) are unknown. But by applying the method described above, it is possible to restore the values of parameters $\alpha$ and ${\lambda }_0$ in the model (5) on the basis of experimental data on RVA variations in the accumulation chamber of the INSR point, and then to estimate the values of $q^{\left(1\right)}$ before the earthquake and $q^{\left(2\right)}$ after the earthquake by (11).

Remark 13.

Algorithms realizing solutions of direct problems: by ODE-model (2) and hereditary α-model (5) underlying the IP-LB algorithm, make calculations in <<meter/hour>> values. This is justified by the fact that the actual volume of the storage chamber (∼0.01 [m${}^3$]) and the frequency of RVA registration (2 cycle./hour). But for convenience of perception, the other parameters of the models (2) and (5) were recorded according to the international system of units SI. Therefore, further in (Figure 5 and Figure 6) parameter values $T,t,h$ in [h], and the recovered AER value $X_1={\lambda }_0$ in [h${}^{-1}$]. Hence, to estimate the RFD (11), we need to divide ${\lambda }_0$ by 3600 to go to the scale [s${}^{-1}$].

Then for the camera in (Figure 2) RFD [Bq/$m^2$s] will be calculated by the formula:

$$q={\displaystyle \frac{A_{max}{\lambda }_{0}V}{3600\Pi }},$$

(12)

7. Modeling Results

Next, we will present the results of IP-LB solution to restore the optimal values of $X_0=\alpha$, $X_1={\lambda }_0$ of the hereditary $\alpha$-model RVA (5) from known experimental RVA data.

Remark 14.

The restore of ${\lambda }_0$ is strongly influenced by the choice of its control parameters in the IP-LB algorithm, i.e., $X_1^{\left(0\right)}$ – the initial specified approximation and ${\delta }_1$ – the initial specified increment of $X_1^{\left(0\right)}$. It was observed that if the initial approximation and its increment are set $\pm 1\%$ of the estimate of the maximum value of the parameter to be recovered, the efficiency and accuracy of IP-LB will be the highest based on the estimates of the similarity coefficients for the described data samples.

The results in (Figure 5 and Figure 6) obtained with the following values of control parameters: $\left[X_0^{\left(0\right)},{\delta }_0\right]=\left[0.05,0.01\right]$ for the index $\alpha$; and $\left[X_1^{\left(0\right)},{\delta }_1\right]=\left[0.0025,0.0005\right]$ for the coefficient ${\lambda }_0$. The values of the control parameters c, v, $\Sigma$ described above for the IP-LB algorithm are shown in the figures.

The degree to which the model curves agree well with the experimental data was evaluated by $\sigma$ – Pearson correlation coefficient [41] with the processed experimental RVA data used to solve inverse problems. The key parameters characterizing the obtained results are reduced to <<meter/second>> according to the international system of units SI and are summarized in the (Table 1).

8. Contribution to research and funding

Data from the KRMR point were provided by Makarov E.O., senior researcher of the laboratory of acoustic and radon monitoring of the Kamchatka branch of the Federal Research Center <<Unified Geophysical Service of the Russian Academy of Sciences>> Petropavlovsk-Kamchatsky, Russia.

The work was supported by the Ministry of Education and Science of Russia (within the framework of the state task No. 075-00682-24) and with the use of data obtained at the unique scientific installation <<Seismoinfrasound complex for monitoring of the Arctic cryolithozone and the complex of continuous seismic monitoring of the Russian Federation, neighbouring territories and the world>>.

All calculations related to the solution of direct and inverse problems by RVA models, as well as calculations on data processing, were performed in the PRPHMM 1.0 software package in MATLAB [47] version R2023b for GNU/Linux Ubuntu Desktop 22.04. The PRPHMM 1. 0 is developed within the framework of the project <<Modeling of dynamic processes in geospheres taking into account heredity>> at the expense of the grant of the Russian Science Foundation № 22-11-00064 (head Parovik R.I.) implemented at the Institute of Cosmophysical Research and Radio Wave Propagation, Far East Branch of the Russian Academy of Sciences, Paratunka village, Russia.

Visualization in this article was performed by Tverdyi D.A. using the software package FEVO 1.0, developed also in the script language Gnuplot 6.0 [48] for GNU/Linux Ubuntu Desktop 22.04. The software package FEVO 1. 0 is developed within the framework of the project ‘Development of a software package for modeling and analysis of volumetric activity of radon as a precursor of strong earthquakes in Kamchatka’ under the grant of the Russian Science Foundation № 23-71-01050 (head Tverdyi D.A.) implemented at the Institute of Cosmophysical Research and Radio Wave Propagation, Far East Branch of the Russian Academy of Sciences, Paratunka village, Russia.

9. Conclusions

The scientific novelty of the presented article consists in the fact that for the first time an hereditary mathematical model has been applied to estimate the change of radon flux density under the change of stress-strain state of the medium in the conditions of a strong earthquake. Moreover, for the first time estimates of radon flux density were given before and after a strong earthquake.

From the results presented in (Table 1) it can be seen that according to the proposed hereditary $\alpha$-model RVA, the radon flux density increased by $12.75\%$, while based on the ODE-model RVA radon flux density estimate shows an increase of $5.39\%$.

The performed estimations also showed that the radon flux density increased after the earthquake, which is probably related to the impact of seismic waves on the ground in the area of installation of the accumulation chamber, as well as with the processes of relaxation of the geo-medium, occurring after the seismic event. This impact led to a change in rock permeability in the area of installation of the accumulation chamber and an increase in the total diffusive-convective flow of geogas. Postseismic anomalous changes in RVA associated with the impact of seismic waves were recorded earlier in Kamchatka after the strong Zhupanovskoe earthquake (30 January 2016, $M_w=7.2$) [3]