Development of a Distributed Control System for the Hydrodynamic Processes of Aquifers, Taking into Account Stochastic Disturbing Factors

Abstract

Uncontrolled, intensive extraction of water from mineral water deposits can lead to negative consequences, the penetration of sewage beyond the water-resistant layers, the destruction of the structure of aquifers and the further loss of the mineral water source. This paper deals with the issues of control of hydrodynamic processes under random impacts in the aquifers of mineral water deposits using mathematical models. The analysis of the technological process of mineral water extraction and statistical analysis of retrospective monitoring data of the deposit were carried out in order to confirm the hypothesis of the stationarity of the processes under consideration. A mathematical model of hydrogeological processes occurring under random disturbances has been constructed. The influence of such disturbances is expressed in deviations of the head in the aquifer from the desired value. Controlling consists of maintaining the head of aquifers at a given level in order to preserve the hydro-mineral base of the region. Lumped and distributed controllers are proposed as corrective blocks, the parameters of which are calculated by frequency methods of controller synthesis. The use of a closed loop control system minimizes the influence of random effects.

1. Introduction

To date, there is active production of heavy water in mining companies, and the demand for it continues to grow. At the same time, uncontrolled intensive exploitation can cause negative changes in the aquifer complex: the penetration of sewage outside the aquifers or the disruption of the structure of aquifers. Therefore, it is necessary to consider hydrogeological processes as control objects with their inherent elements and control links, taking into account random effects occurring in a hydrogeological object. Random factors affecting the study object are: resource consumption per debit of water consumption, stochastic change in water flow rate, increased seismic activity [1,2,3,4,5,6].

When describing a mineral water deposit, its extent in space should be taken into account, as a result of which it is required to represent the control object as a system with distributed parameters. In addition, when studying the system, it is necessary to consider the influence of random factors on hydrogeological processes.

The most detailed basic information about the methods for constructing mathematical models of hydrological processes is found in the works of Gavich I.K., Antonov V.V., Lomakin E.A., Konosavsky P.K. and others. Their work has made a great contribution to the development of methods for constructing two-dimensional models. However, these works in their pure form cannot be used to build three-dimensional models of spatially distributed objects. Methodology for the analysis of distributed objects, in particular representation in the frequency domain, was carried out by such scientists as Butkovsky A.G., Sirazetdinov T.K. and others [7]. Among the studies of recent years in the field of control systems with distributed parameters, a significant role was played by the studies of Malkov A.V. and Pershin I.M., who, in their works solve the problems of improving the theoretical foundations of the management of hydro-geological objects [8]. In the study of stochastic systems, there are the works of many scientists, among which are Ostrem K., Evlanov L.G., Konstantinov V.M. and others. The listed works in the field of controls, both lumped [9,10,11] and distributed systems [12,13,14,15], do not address the analysis of stochastic processes in distributed control objects and specifically the control of hydrodynamic processes, taking into account the random nature of the change in object parameters.

Thus, the main tasks aimed at increasing the economic efficiency of mineral water extraction while minimizing the environmental threats [16,17] are the analysis of hydrogeological processes under random impacts, the search for control methods and algorithms, and the laws of optimal functioning of the technological process.

2. Topical Area Analysis and Formulation of the Problem

Analyzing hydrolithospheric processes as a control object, it is possible to reveal the structure of the interaction of aquifers in a hydrogeological object and identify parameters that may have a random change in time and complicate the process of maintaining pressure in the aquifer at a given level. Allowing the head to drop below the maximum acceptable rate, water from the upper aquifers begins to abundantly seep into the lower ones, especially from the groundwater layer. The quality of mineral water in the exploited aquifer worsens. This will lead to a high risk of losing the exploited source of mineral water entirely.

The exploited area of the Kislovodsk mineral water deposit is considered as the object of study. Monitoring in observation wells showed systematic fluctuations in water pressure over time. There can be several reasons for this: damage to equipment, precipitation, changes in porosity, and human factors.

These systematic fluctuations due to various factors can be described and considered as uncertain or random processes, and these processes should be included when designing a mathematical model of hydrogeological objects. Additionally, it is worth considering seismic activity, which also has an uncertain or random impact. The assumed random factors have a stochastic effect on the flow rate of water in the aquifer and on the level of head drawdown. Analyzing the dynamics of levels and their systematic fluctuations, it is possible to suggest some limitation for the mathematical model of the object, that the simulated random processes will be considered as stationary. To test this hypothesis a statistical analysis of piezometric level monitoring data in observation wells of the considered area of the field was performed.

The purpose of this study was to develop a control system for the hydrolithospheric process under random impacts, aimed at maintaining the quality and reserves of the mineral water of the deposit by controlling the piezometric level in the exploited aquifer.

It is required to carry out a detailed analysis of a hydrogeological object with an analysis of its inherent random processes that affect the hydrodynamic parameters during the operation in the field. This requires system-analytical modeling of the technical and physical elements of the object of study at different levels of management with the further development of a mathematical model that takes into account random disturbances in the geofiltration processes of aquifers of a mineral water deposit. The mathematical model is supplemented with initial and boundary conditions and verified with data obtained experimentally in a real mineral water production facility.

3. Study Area Description

The study area is located in the territory of two regional taxonomic units of hydrogeological zoning of the 1st order (according to L.K. Ostrovsky): the Crimean Ciscaucasian province of complex basins of reservoirs and block-stratal waters and the Crimean Ciscaucasian province of complex basins of reservoir-block, reservoir, karst-block, vein-block and cover-flow groundwater. The boundaries between them pass along the system of Ciscaucasian deep faults in the basement. Within the first province, two basins of groundwater of the second order are located: there are the Azov–Kuban and East Ciscaucasian. And within the second province there are the Greater Caucasus and the Central Caucasian.

The territory of the Greater Region of the Caucasian Mineral Waters (CMW) is located in the junction zone of the indicated hydrogeological provinces along the boundaries of the Greater Caucasus and East Ciscaucasian artesian basins. The Kislovodsk mineral water deposit is located within the Bolshekavkazsky artesian basin.

In geological and structural terms, the territory of the Greater CMW region is part of the northern wing of the Greater Caucasus meganticlinorium, known as the North Caucasian monocline, and is located in the region of the meridional Mineralovodsky uplift. This uplift has the appearance of a wide and gently sloping anticline, gradually plunging to the north. The arch of the uplift, up to 60 km wide, is complicated in places by transverse folds and faults in the longitudinal and transverse directions. The eastern and western wings of the Mineralnye Vody uplift plunge into the Terek–Kuma and Kuban depressions, respectively.

In tectonic terms the study area is a part of the Alpine folding, within which two structural stages are distinguished: the crystalline basement and the sedimentary cover.

The crystalline basement of the Lower Paleozoic age is represented by dislocated crystalline greenish–gray chlorite–albite and mica-chlorite schists. The basement is divided by a system of tectonic discontinuities of northeastern strike with an amplitude of faults and reverse faults up to 20–50 m, and is complicated by felsic intrusions, mainly pink coarse-grained granites.

The general patterns of occurrence of aquifers are clear from the geological and hydrogeological profile (Figure 1).

The oldest rocks on the territory of the Kislovodsk deposit are Lower Paleozoic (PZ₁) crystalline schists and Upper Paleozoic (PZ) granites. They were discovered by wells at depths ranging from 80 m in the southern and up to 400 m in the northern parts of the deposit.

Jurassic system (J): The sediments of the Jurassic system lie with angular and stratigraphic unconformity on the Paleozoic sediments and are represented by the Toarcian stage of the Lower Jurassic and the Tithonian stage of the Upper Jurassic. The remaining stratigraphic sediments of the Jurassic system are eroded.

Toarcian Stage (J₁t): Toarcian sediments are found only in the upper southern part of the Berezovskaya Gully. To the north, they are cut off by the pre-Tithonian transgression.

This provides a basis for believing that the Toarcian sediments, lying on the surface of the crystalline basement and performing irregularities in its relief, formed lenses during subsequent erosion, one of which lies in the area of the Berezovsky section of the Kislovodsk deposit. The sediments are represented by mudstones and sandstones which are up to 70 m thick.

Tithonian Stage (J₃tt): The sediments of the Tithonian stage hide on the eroded surface of the basement. They were discovered by wells at depths from 100 to 400 m, increasing in a northerly direction. The Tithonian variegated sequence consists of loose sandstones, clays, and granite gruss. The sediments of the Tithonian stage are divided by researchers into three sequences: the lower one, which is dominated by coarse-grained variegated sandstones with rare interlayers of gruss and clays; the middle one, consisting of granite gruss with interlayers of sandstones and clays, and the upper one, composed of red–brown clayey sandstones with gypsum inclusions. The thickness of the sediments of the Tithonian stage is inversely related to the relief of the top of the Paleozoic basement. In areas where the basement surface forms depressions, the thickness of the variegated stratum increases, and, conversely, in the elevated areas of the basement, it is minimal.

Cretaceous system (K): The sediments of the Cretaceous system are represented by the lower section. The Upper Cretaceous section in the area of the site is completely eroded. The Lower Cretaceous sediments are represented by the Valanginian, Hauterivian, Barremian, Aptian, and Albian stages.

Valanginian Stage (K₁v). Valanginian sediments, represented in the section by dolomitic limestones, overlie the variegated Tithonian sequence. They are exposed on the surface in a small area in the southern part of the deposit, along the valleys of the rivers Berezovaya, Kabardinka and Olkhovka.

Dolomitized limestones were discovered in the deposit at depths from 5 to 300 m, increasing from south to north. Their thickness is quite consistent over the area and varies within insignificant limits—from 80 to 90 m.

Hauterivian Stage (K₁g): The sediments of the Hauterivian stage lie on the dolomitic limestones of the Valanginian and are important in terms of sandy-clay composition. The sediments were formed by a number of shell limestones and argillaceous sandstones with limestone interbeds and outcrops in the main and northern parts of the field, as well as in rarely elevated areas.

The total thickness of the Hauterivian sediments is about 50 m, and like the underlying sediments, have a slight slope to the northeast.

Barremian (K₁br): The sediments of the Barremian Stage are represented by a sequence of dense oolitic limestones interbedded with limestone–argillaceous sandstones. Above are calcareous sandstones with a small number of thin layers of oolitic limestones. They are overlain by argillaceous siltstones, which in some places change into strongly argillaceous dark brown sandstones and higher up are replaced by grayish-yellow fine-grained, argillaceous sandstones. The section of the Barremian sediments ends with red ferruginous sandstones. Sediment thickness about 40 m.

Aptian Stage (K₁a): The Aptian stage is represented by clays (lower substage) overlain by sands and sandstones (upper substage). The total thickness of the sediments of the Aptian stage is 245 m.

Albian Stage (K₁al): The sediments of the Albian stage are found in the north and west of the described territory. In the lower part, they are represented by sands and sandstones (lower substage), and in the upper part, by black calcareous clays (upper substage). The thickness of the Albian Stage is 150–200 m.

To the north, the Lower Cretaceous sediments subside under the Upper Cretaceous sediments.

Cenozoic era (KZ): In the area of the Kislovodsk deposit, it is represented by eluvial, deluvial and alluvial formations of the Upper Quaternary system (Q_IV), which are distributed mainly along the thalwegs of permanent and temporary watercourses. Eluvial sediments are widespread in the watershed areas north and south of the Berezovaya, Alikonovka, Podkumok rivers. They are represented by yellowish-gray compacted loams, containing crushed sandstones mainly in the lower part of the section.

Deluvial sediments are most developed along the slopes of the valleys of the rivers mentioned above, especially in the lower part of the slopes. In the upper part of the slopes, deluvial sediments are insignificantly thick, and often completely absent. They are represented mainly by yellowish-gray and reddish-brown loams with a large number of fragments of dolomitic limestones of different sizes: from small gravel to huge blocks. There are especially many blocks in the middle part of the valleys, where their slopes are completely composed of dolomitic limestones. Here, these blocks often reach the riverbeds and, cluttering them up, change the direction of the river flow. The thickness of deluvial sediments varies from 0.2 to 50 m.

Alluvial sediments fill the bottom of river valleys and ravines flowing into them. Here they are represented by alternating clays, loams, sand and poorly rounded pebbles and gravel. Boulders of various sizes are also found in alluvial sediments. In the region of the river Berezovaya at the foot of the right slope there are often silt-loam sediments, as well as buried soils and travertines. The mountainous nature of the river channels and their destructive activity during flood periods create unfavorable conditions for the accumulation of precipitation, which generally determines the insignificant thickness of alluvial sediments—no more than 5–10 m, with the exception of some areas where it reaches 40 m.

In addition to the geological characteristics of the deposit, information about the mineral water reserves of the research aquifer is in

Table 1. Information about the mineral water reserves of the Kislovodsk deposit.

Well №	Exploitation Reserves, m3/Days	Drilling Year
23	3	1955
8	2	1963

. There are some producing wells, which are produced from water of the Tithonian aquifer.

The state of the aquifer exploitation and piezometric monitoring are shown in Figure 2.

4. Research Methods

Analysis of modern research has made it possible to determine a set of effective methods for modeling hydrogeological objects, which are most often used to solve problems in this subject area [8,13]. As a result of the presented method’s effectiveness, the method using a mathematical description of the hydrogeological object in the form of a three-dimensional geo-filtration model was chosen for analysis of the further research. The experiment carried out on the prototype of the hydrogeological object model showed convergence values of the real data with the model ones by 80 to 90%. Errors in the modeling of geo-filtration processes can lead to errors in controlling the aquifer system’s exploitation. Due to its effectiveness, this method of modeling is used for synthesis distributed control systems.

Before proceeding to the mathematical model of the object, studies of the piezometric level in observation wells should be conducted. It is necessary in order to put forward a hypothesis about the stationarity of the considered random processes and take these processes into account when modeling the object. Most of the methods of analysis and forecasting of dynamic processes are applicable to stationary time series. After analyzing the hydrogeological object of study, it was concluded that a sufficient method for studying stationarity would be to evaluate the series relative to its central trend. Often, to detect the non-stationarity of a series, it is enough to consider the process plot. At the same time, for statistically valid proof of stationarity, it is necessary to use special methods, such as compiling a correlogram of the process and splitting into successive groups of observations to compare them with each other.

For the analysis of a series of dynamics of the piezometric level from observation well No. 87 in the field using the correlogram, it is necessary to supplement the initial series of dynamics with series from the levels of this series shifted in time (

Table 2. The initial series of level dynamics (in meters) with shifted series for statistical processing.

Yt	Yt-1	Yt-2	Yt-3	Yt-4	Yt-5	Yt-6	Yt-7
821.05
821.12	821.05
821.29	821.12	821.05
820.96	821.29	821.12	821.05
820.89	820.96	821.29	821.12	821.05
820.89	820.89	820.96	821.29	821.12	821.05
820.89	820.89	820.89	820.96	821.29	821.12	821.05
820.75	820.89	820.89	820.89	820.96	821.29	821.12	821.05
820.75	820.75	820.89	820.89	820.89	820.96	821.29	821.12
820.75	820.75	820.75	820.89	820.89	820.89	820.96	821.29
821.22	820.75	820.75	820.75	820.89	820.89	820.89	820.96
821.67	821.22	820.75	820.75	820.75	820.89	820.89	820.89
822.05	821.67	821.22	820.75	820.75	820.75	820.89	820.89
822.44	822.05	821.67	821.22	820.75	820.75	820.75	820.89
822.08	822.44	822.05	821.67	821.22	820.75	820.75	820.75
821.55	822.08	822.44	822.05	821.67	821.22	820.75	820.75
821.58	821.55	822.08	822.44	822.05	821.67	821.22	820.75
822.22	821.58	821.55	822.08	822.44	822.05	821.67	821.22
822.91	822.22	821.58	821.55	822.08	822.44	822.05	821.67
823.07	822.91	822.22	821.58	821.55	822.08	822.44	822.05
823.18	823.07	822.91	822.22	821.58	821.55	822.08	822.44
823.52	823.18	823.07	822.91	822.22	821.58	821.55	822.08
823.94	823.52	823.18	823.07	822.91	822.22	821.58	821.55
824.06	823.94	823.52	823.18	823.07	822.91	822.22	821.58
824.06	824.06	823.94	823.52	823.18	823.07	822.91	822.22
823.85	824.06	824.06	823.94	823.52	823.18	823.07	822.91
822.86	823.85	824.06	824.06	823.94	823.52	823.18	823.07
823.2	822.86	823.85	824.06	824.06	823.94	823.52	823.18
822.97	823.2	822.86	823.85	824.06	824.06	823.94	823.52
822.25	822.97	823.2	822.86	823.85	824.06	824.06	823.94
822.15	822.25	822.97	823.2	822.86	823.85	824.06	824.06
821.79	822.15	822.25	822.97	823.2	822.86	823.85	824.06
821.91	821.79	822.15	822.25	822.97	823.2	822.86	823.85
822.25	821.91	821.79	822.15	822.25	822.97	823.2	822.86
822.34	822.25	821.91	821.79	822.15	822.25	822.97	823.2
822.29	822.34	822.25	821.91	821.79	822.15	822.25	822.97

). Calculating the correlation coefficients of the initial series Y_t with other ones, a resulting correlogram can be presented in the following form in Figure 3. Since the first-order autocorrelation coefficient turned out to be the highest, the time series contains a trend. As the delay increases, the correlation monotonically decreases, which indicates the stationarity of the series, the probabilistic properties of which do not change with time. To confirm the results, the second method was used for estimating the stationarity of the series.

To begin, data must be selected and grouped. Thus, for the observation period from 2008 to 2020 (150 observations) we divided the series into three equal segments (50 values each). Differences in the means of the series are insignificant (t = 0.96, df = 98, p > 0.05), where t is Student’s t-test, df is the degree of freedom, and p is the level of significance. The data was grouped by 100 observations; as an example two segments were chosen (2004–2012 and 2012–2020). The differences in the means of the series are also insignificant (t = 0.64, df = 198, p > 0.05). Thus, the stationarity of the time series according to the observation well data is proven. A number of other control wells of the field were analyzed using the same methods.

The conducted studies were systematized and presented in the form of a methodology that was recommended to be used for the statistical processing of data on observation wells of a hydrogeological object in order to assess the stationarity of processes. Having conducted a study of the effect of random disturbances in the reservoir, we moved on to a system analysis of a hydrogeological object for the subsequent development of a mathematical model.

To simplify the description and to highlight the disturbing and input effects, a visual model is proposed in the form of a block diagram (Figure 4), which displays not only the effects on the exploited aquifer, but also the interaction between the layers.

In many literary sources [8,18,19,20,21], the basic principles of mathematical modeling of aquifer systems and various modeling options are given. Preference was given to describing the object using a three-dimensional model. This type of mathematical model describes filtration processes in aquifers and interactions between them through various coefficients. The mathematical model of the central area of the Kislovodsk deposit with four producing wells and four observation wells can be written as a system of differential equations that determine the geofiltration processes of the aquifer system. Let ${\mathrm{h}}_1$ be the head in the groundwater and ${\mathrm{H}}_{\mathrm{i}}$ be the head in i-th aquifer (i = 2,3,4), then the mathematical model of the field can be written as follows:

$$\frac{\partial {\mathrm{h}}_1(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathsf{\tau }}=\frac{1}{{\mathsf{\eta }}_1}\left({\mathrm{k}}_{1,\mathrm{x}}\frac{{\partial }^2{\mathrm{h}}_1(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{x}}^2}+{\mathrm{k}}_{1,\mathrm{y}}\frac{{\partial }^2{\mathrm{h}}_1(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{y}}^2}+{\mathrm{k}}_{1,\mathrm{z}}\frac{{\partial }^2{\mathrm{h}}_1(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{z}}_1{}^2}\right);$$

(1)

$$\begin{array}{c}\frac{\partial {\mathrm{H}}_2(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathsf{\tau }}=\\ =\frac{1}{{\mathsf{\eta }}_2}\left({\mathrm{k}}_{2,\mathrm{x}}\frac{{\partial }^2{\mathrm{H}}_2(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{x}}^2}+{\mathrm{k}}_{2,\mathrm{y}}\frac{{\partial }^2{\mathrm{H}}_2(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{y}}^2}+{\mathrm{k}}_{2,\mathrm{z}}\frac{{\partial }^2{\mathrm{H}}_2(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{z}}_2{}^2}\right)-\\ -({\overline{\mathrm{F}}}_{2,\mathrm{x}}+{\tilde{\mathrm{F}}}_{2,\mathrm{x}})\cdot \frac{\partial {\mathrm{H}}_2(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathrm{x}};\end{array}$$

(2)

$$\begin{array}{c}\frac{\partial {\mathrm{H}}_3(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathsf{\tau }}=\\ =\frac{1}{{\mathsf{\eta }}_3}\left({\mathrm{k}}_{3,\mathrm{x}}\frac{{\partial }^2{\mathrm{H}}_3(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{x}}^2}+{\mathrm{k}}_{3,\mathrm{y}}\frac{{\partial }^2{\mathrm{H}}_3(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{y}}^2}+{\mathrm{k}}_{3,\mathrm{z}}\frac{{\partial }^2{\mathrm{H}}_3(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{z}}_3{}^2}\right)-\\ -({\overline{\mathrm{F}}}_{3,\mathrm{x}}+{\tilde{\mathrm{F}}}_{3,\mathrm{x}})\cdot \frac{\partial {\mathrm{H}}_3(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathrm{x}};\end{array}$$

(3)

$$\begin{array}{c}\frac{\partial {\mathrm{H}}_4(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathsf{\tau }}=\\ =\frac{1}{{\mathsf{\eta }}_4}\left({\mathrm{k}}_{4,\mathrm{x}}\frac{{\partial }^2{\mathrm{H}}_4(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{x}}^2}+{\mathrm{k}}_{4,\mathrm{y}}\frac{{\partial }^2{\mathrm{H}}_4(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{y}}^2}+{\mathrm{k}}_{4,\mathrm{z}}\frac{{\partial }^2{\mathrm{H}}_4(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial {\mathrm{z}}_4{}^2}\right)-\\ -({\overline{\mathrm{F}}}_{4,\mathrm{x}}+{\tilde{\mathrm{F}}}_{4,\mathrm{x}})\cdot \frac{\partial {\mathrm{H}}_4(\mathrm{x},\mathrm{y},\mathrm{z},\mathsf{\tau })}{\partial \mathrm{x}}+(\overline{\mathrm{V}}({\mathrm{y}}_{\mathrm{j}},\mathsf{\tau })+\tilde{\mathrm{V}}({\mathrm{y}}_{\mathrm{j}},\mathsf{\tau }))\cdot \mathsf{\delta }({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{z}}_{\mathrm{j}});\\ 0<\mathrm{x}<{\mathrm{L}}_{\mathrm{x}},0<\mathrm{y}<{\mathrm{L}}_{\mathrm{y}},0<\mathrm{z}<{\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}}\end{array}$$

(4)

where: ${\mathsf{\eta }}_{\mathrm{i}}$—elastic capacity parameter of the i-th aquifer; ${\overline{\mathrm{F}}}_{\mathrm{i},\mathrm{x}}+{\tilde{\mathrm{F}}}_{\mathrm{i},\mathrm{x}}={\mathrm{F}}_{\mathrm{i},\mathrm{x}}$—flow speed of the i-th aquifer, represented as the sum of the deterministic and random components, respectively; $\overline{\mathrm{V}}({\mathrm{y}}_{\mathrm{j}},\mathsf{\tau })+\tilde{\mathrm{V}}({\mathrm{y}}_{\mathrm{j}},\mathsf{\tau })=\mathrm{V}({\mathrm{y}}_{\mathrm{j}},\mathsf{\tau })$—head reduction in producing wells, represented as the sum of the deterministic and random components, respectively; k_i,x, k_i,y, k_i,z—filtration coefficients; ${\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{z}}_{\mathrm{j}}$—coordinates of the location of producing wells (j = 1,2,3,4); $\mathsf{\delta }({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{z}}_{\mathrm{j}})$—parameter that is equal to one when x = x_j, y = y_j, z = z_j, otherwise it is equal to zero.

In the considered mathematical model, the filtration parameter k (m/day) is the most important, since it determines not only the nature of water overflow in the horizon, but also makes it possible to quantify its volume. The indicator of this parameter is determined by the hydro-chemical components of water and the porosity of the studied rock. This parameter is determined using experimental filtration works. Based on measurements at certain points (location of observation wells), the average value of the parameter in the study area is found.

The geo-filtration model is complicated by the fact that it takes into account the random nature of external factors in the form of 10% stochastic noise head reduction at the locations of production wells and at flow speed in the aquifers.

The mathematical model is also supplemented with boundary conditions [22].

Boundary conditions between aquifers (Darcy conditions):

$$\begin{gathered} {\mathrm{h}}_1(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_1},\mathsf{\tau })={\mathrm{h}}_1(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_1},\mathsf{\tau })+{\mathrm{b}}_1\cdot \left({\mathrm{H}}_2(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{h}}_1(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_1},\mathsf{\tau })\right); \\ {\mathrm{H}}_2(\mathrm{x},\mathrm{y},0,\mathsf{\tau })={\mathrm{H}}_2(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{b}}_1\cdot \left({\mathrm{H}}_2(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{h}}_1(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_1},\mathsf{\tau })\right); \\ {\mathrm{H}}_{\mathrm{i}}(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}},\mathsf{\tau })={\mathrm{H}}_{\mathrm{i}}(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}},\mathsf{\tau })+{\mathrm{b}}_{\mathrm{i}}\cdot \left({\mathrm{H}}_{\mathrm{i}}(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{H}}_{\mathrm{i}}(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}},\mathsf{\tau })\right); \\ {\mathrm{H}}_{\mathrm{i}+1}(\mathrm{x},\mathrm{y},0,\mathsf{\tau })={\mathrm{H}}_{\mathrm{i}+1}(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{b}}_{\mathrm{i}}\cdot \left({\mathrm{H}}_{\mathrm{i}+1}(\mathrm{x},\mathrm{y},0,\mathsf{\tau })-{\mathrm{H}}_{\mathrm{i}}(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}},\mathsf{\tau })\right); \end{gathered}$$

(5)

where ${\mathrm{b}}_{\mathrm{i}}$—overflow parameter (i = 1…3).

Lower boundary of the lower aquifer:

$$\partial {\mathrm{H}}_4(\mathrm{x},\mathrm{y},{\mathrm{L}}_{{\mathrm{z}}_4},\mathsf{\tau })/\partial \mathrm{z}=0.$$

(6)

An assumption in the model is that the thickness of the reservoirs along the y coordinate is so high that the impacts from production wells do not significantly affect the reservoir pressure at the horizon boundaries. Then the boundary conditions of the side boundaries are written in the following order:

$$\begin{gathered} {\mathrm{h}}_1(0,\mathrm{y},\mathrm{z},\mathsf{\tau })={\mathrm{h}}_{10}; \\ {\mathrm{h}}_1(\mathrm{x},0,\mathrm{z},\mathsf{\tau })={\mathrm{h}}_1(\mathrm{x},{\mathrm{L}}_{\mathrm{y}},\mathrm{z},\mathsf{\tau })={\mathrm{h}}_{1,0}; \\ {\mathrm{H}}_{\mathrm{i}}(0,\mathrm{y},\mathrm{z},\mathsf{\tau })={\mathrm{H}}_{\mathrm{i}0}; \\ \partial {\mathrm{H}}_{\mathrm{i}}({\mathrm{L}}_{\mathrm{x}},\mathrm{y},\mathrm{z},\mathsf{\tau })/\partial \mathrm{x}=0; \\ {\mathrm{H}}_{\mathrm{i}}(\mathrm{x},0,\mathrm{z},\mathsf{\tau })={\mathrm{H}}_{\mathrm{i}}(\mathrm{x},{\mathrm{L}}_{\mathrm{y}},\mathrm{z},\mathsf{\tau })={\mathrm{H}}_{\mathrm{i}0}; \end{gathered}$$

(7)

where ${\mathrm{h}}_{10}$, ${\mathrm{H}}_{\mathrm{i}0}$—initial states of unperturbed aquifers.

The geometric and physical data of the deposit (in meters) are given in [22]. To assess the dynamics of hydrodynamic processes, using a difference scheme, the object model is converted to a discrete form, since for the system of equations of geo-filtration processes in this setting taking into account the given boundary and initial conditions of the modeled area, it is difficult to obtain an analytical solution [23,24,25,26,27].

To compile a discrete form of the above-described mathematical model of the field, the area under study was divided by some grids with the following parameters: $\Delta \mathrm{x}={\mathrm{L}}_{\mathrm{x}}/({\mathrm{N}}_{\mathrm{x}}-1)$, $\Delta \mathrm{y}={\mathrm{L}}_{\mathrm{y}}/({\mathrm{N}}_{\mathrm{y}}-1)$, $\Delta {\mathrm{z}}_{\mathrm{i}}={\mathrm{L}}_{{\mathrm{z}}_{\mathrm{i}}}/({\mathrm{N}}_{{\mathrm{z}}_{\mathrm{i}}}-1)$—discrete steps to the respective coordinate axes (i = 1…4), where ${\mathrm{N}}_{\mathrm{x}}$, ${\mathrm{N}}_{\mathrm{y}}$, ${\mathrm{N}}_{{\mathrm{z}}_{\mathrm{i}}}$—the number of discrete points (in our case ${\mathrm{N}}_{\mathrm{x}}=8$, ${\mathrm{N}}_{\mathrm{y}}=8$, ${\mathrm{N}}_{{\mathrm{z}}_{\mathrm{i}}}=7$).

For the exploited lower aquifer, the equation in discrete form will be:

$$\begin{array}{c}\frac{\Delta {\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }}}{\Delta \mathsf{\tau }}=\frac{1}{{\mathsf{\eta }}_4}({\mathrm{k}}_{3,\mathrm{x}}\frac{{\mathrm{H}}_{4,\mathsf{\eta }-1,\mathsf{\gamma },\mathsf{\xi }}-2\cdot {\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }}+{\mathrm{H}}_{4,\mathsf{\eta }+1,\mathsf{\gamma },\mathsf{\xi }}}{\Delta {\mathrm{x}}^2}+{\mathrm{k}}_{4,\mathrm{y}}\frac{{\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma }-1,\mathsf{\xi }}-2\cdot {\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }}+{\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma }+1,\mathsf{\xi }}}{\Delta {\mathrm{y}}^2}+\\ +{\mathrm{k}}_{4,\mathrm{z}}\frac{{\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }-1}-2\cdot {\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }}+{\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }+1}}{\Delta {\mathrm{z}}_4{}^2})-({\overline{\mathrm{F}}}_{4,\mathsf{\eta }}+{\tilde{\mathrm{F}}}_{4,\mathsf{\eta }})\cdot \frac{{\mathrm{H}}_{4,\mathsf{\eta },\mathsf{\gamma },\mathsf{\xi }}-{\mathrm{H}}_{4,\mathsf{\eta }-1,\mathsf{\gamma },\mathsf{\xi }}}{\Delta \mathrm{x}}+({\overline{\mathrm{V}}}_{\mathsf{\gamma }}+{\tilde{\mathrm{V}}}_{\mathsf{\gamma }})\cdot {\mathsf{\delta }}_{\mathsf{\eta }0,\mathsf{\gamma }0,\mathsf{\xi }0};\\ 2<\mathsf{\eta }<{\mathrm{N}}_{\mathrm{x}}-1;2<\mathsf{\gamma }<{\mathrm{N}}_{\mathrm{y}}-1;2<\mathsf{\xi }<{\mathrm{N}}_{\mathrm{z}4}-1.\end{array}$$

(8)

The remaining equations are written in a similar way.

Using a discrete model, software was developed that numerically simulates the hydrodynamic processes of aquifers under random impacts. The simulation results are used to analyze the dynamic characteristics of the object and the subsequent synthesis of the control system for the mineral water production process.

The graph of the object’s response to the first spatial harmonic of the input action, obtained at the location of the second observation well, and the dynamics of the head deviation under the effects of random disturbances are illustrated in Figure 5.

The head deviation in the steady state reaches 0.12 m, which can lead to the penetration of sewage into the exploited aquifer. The problem is to control the hydrodynamic process and reduce the effect of random disturbances on maintaining the head of the exploited reservoir at a given level in order to preserve the reserves and quality of mineral water.

The frequency method for the synthesis of distributed controllers was developed for distributed objects, the models of which are determined using experimental data [8]. This method has shown itself effective in controlling various objects of thermal and electromagnetic natures. Additionally, this method can be applied to controlling hydrogeological research objects, the model of which is obtained using experimental data.

The structure of a distributed control system can be presented as in Figure 6.

The control problem is solved by two methods: by the synthesis of lumped controllers from the position of considering the object as a multidimensional system and by the synthesis of distributed controllers using frequency synthesis methods of systems with distributed parameters. Using approximating functions that describe the reactions of the object to the input action [28,29,30,31,32], the following controllers are synthesized.

Transfer function of calculated PID controller:

$$\mathrm{R}(\mathrm{s})={\mathrm{K}}_{\mathrm{p}}+\frac{1}{{\mathrm{T}}_{\mathrm{i}\mathrm{n}}\cdot \mathrm{s}}+{\mathrm{T}}_{\mathrm{d}}\cdot \mathrm{s}=0.413+\frac{1}{158437.673\cdot \mathrm{s}}+0.0002\cdot \mathrm{s}$$

(9)

where K_p—proportional gain, T_d—derivative time constant, T_in—integral time constant.

The transfer function of the calculated distributed controller:

Transfer function of calculated PID controller:

$$\begin{gathered} \mathrm{W}\left(\mathrm{y},\mathrm{s}\right)={\mathrm{E}}_1\cdot \left[\frac{{\mathrm{n}}_1-1}{{\mathrm{n}}_1}-\frac{1}{{\mathrm{n}}_1}{\nabla }^2\right]+{\mathrm{E}}_4\cdot [\frac{{\mathrm{n}}_4-1}{{\mathrm{n}}_4}-\frac{1}{{\mathrm{n}}_4}{\nabla }^2]\cdot \frac{1}{\mathrm{s}}+{\mathrm{E}}_2\cdot [\frac{{\mathrm{n}}_2-1}{{\mathrm{n}}_2}-\frac{1}{{\mathrm{n}}_2}{\nabla }^2]\cdot \mathrm{s}; \\ \mathrm{W}\left(\mathrm{y},\mathrm{s}\right)=25.139\cdot \left[\frac{0.001}{1.001}-\frac{1}{1.001}{\nabla }^2\right]+0.0001\cdot \left[\frac{0.015}{1.015}-\frac{1}{1.015}\cdot {\nabla }^2\right]\frac{1}{\mathrm{s}}+1079.880\cdot \mathrm{s}, \end{gathered}$$

(10)

where ${\mathrm{E}}_{\mathrm{i}}$,${\mathrm{n}}_{\mathrm{i}}$—gains and weight coefficients, respectively.

5. Results

By forming control actions and representing them in a discrete form, algorithms for controlling hydrodynamic processes were obtained. This algorithm was implemented in software that simulates the behavior of a closed loop control system. Figure 7 and Figure 8 show the transient processes of the system using the developed controllers.

The results of simulation of closed loop systems show that the parameters of the synthesized controllers are found correctly, and the systems ensure that the pressure is maintained at a given level. But it is worth highlighting the advantage of a distributed controller over a lumped one, since the overshoot is much less (σ_dis = 20% < σ_pid = 53%), and the system reaches a steady state faster.

It is worth noting how small the effects of random disturbances on the control system are. Comparison of the random signal at the output of the system and at the output of the object before applying the distributed controller in Figure 9 illustrates that the random effect was reduced (from ±0.15 m to ±0.005 m).

Geofiltration parameters of the simulated area of the field may change over time. Therefore, there was a check of the stability of a distributed control system to parametric disturbances. The developed control system did not lose its operability with a decrease or increase in filtration coefficients by ten times.

These methods can be applied in the regulation of water resources considering the effect of random disturbing factors [15,16]. In addition, the obtained results can be also used in related areas, where objects can be represented in the context of a class of control systems with distributed parameters [33,34,35,36,37,38,39,40,41,42], for human–machine systems [43,44,45,46], and especially for the oil and gas industry [47,48,49,50,51]. The use of these methods will expand the possible ways of solving some economic issues [52,53,54].

6. Discussion

Recently, there has been an increase in the extraction of mineral waters for industrial bottling and for resort purposes. Irrational water intake can lead to a deterioration in the quality of produced water. From the experimental data of the producing company, there are systematic fluctuations in water head. There can be several reasons: damage to equipment, precipitation, changes in porosity, and human factors. This deposit is located near Mount Elbrus, so seismic influence is not excluded. In addition, anthropogenic impacts play an important role. Regarding the latter, the following should be added: since 2010, there has been active development, both in the central part of the city Kislovodsk and its outskirts. Due to the lack of central sewerage, many people use cesspools, which, under the influence of precipitation, leaks into groundwater. Since 2012, contamination of mineral springs has been systematically observed. With a sufficiently high reservoir pressure, contamination is minimized. But due to uncontrolled water withdrawal, a drop in pressure (head) occurs, which leads to an aggravation of the situation.

All these disturbances in this study are considered as a systematic (or possible) effect on the flow rate of water in the aquifer and the head in the wells. As a result, the piezometric level deviates from the desired value. If the deviations are large, then the risk of flow of contaminants into the exploited aquifer increases. This problem can be solved by reducing water production and increasing the set head value, but at the same time, the profit of the producing company will also decrease. However, the use of a synthesized distributed control system shows minimization of the influence of random processes, thus increasing the profitability.

The advantage of the method of considering random effects in the mathematical model makes it possible to reduce their influence. However, fluctuations can be intermittent and take various forms. Therefore, a feedback control system was applied which is able to control the process under uncertain changes. For example, several experiments were carried out with a sharp increase in the parameters of the model. With such uncertain influences, the system is stable and goes into a static mode, but with a longer regulation time.

The advantage of distributed controller, as was mentioned, is less overshoot than with a lumped controller. When using a lumped controller, the head drops significantly lower than when using a distributed controller. This fact increases the risk of flow of contaminants into the aquifer from the upper ones (due to the head difference). Speed plays an important role, because it shows how quickly the control system stabilizes the production process under uncertain disturbances.

However, the advantage of lumped controllers is a simpler implementation compared to distributed controllers.

7. Conclusions

The paper considers methods for the study of hydrogeological objects, distributed control systems and systems under random influences. On the basis of the study performed, conclusions were drawn about the effectiveness of the use of controllers to regulate the processes of pumping mineral water from the exploited aquifer. The main tasks implemented in the framework of the study included:

• A mathematical model that was developed for the study area of the Kislovodsk mineral water deposit, which took into account the impact of random factors on the flow speed and water debits in aquifers.
• A closed loop system was developed for controlling the hydrolithospheric process under random influences. In addition, the use of distributed controllers minimized the deviation of a random variable at the output of the control system.
• Software was developed for modeling a distributed control system of the hydrodynamic process in aquifers under random impacts.