Do ARMA Models Provide Better Gap Filling in Time Series of Soil Temperature and Soil Moisture? The Case of Arable Land in the Kulunda Steppe, Russia

: The adoption of climate-smart agriculture requires the comprehensive development of environmental monitoring tools, including online observation of climate and soil settings. They are often designed to measure soil properties automatically at different depths at hour or minute intervals. It is essential to have a complete dataset to use statistical models for the prediction of soil properties and to make short-term decisions regarding soil tillage operations and irrigation during a vegetation period. This is also important in applied hydrological studies. Nevertheless, the time series of soil hydrological measurements often have data gaps for different reasons. The study focused on solving a problem of gap-ﬁlling in hourly time series of soil temperature and moisture, measured at the 30 cm depth using a weighted gravitation lysimeter station while meteorological data were recorded simultaneously by a weather station. The equipment was installed in the Kulunda Steppe in the Altai Krai, Russia. Considering that climate conditions affect soil temperature and moisture content directly, we did a comparative analysis of the gap-ﬁlling performance using the three imputation methods—linear interpolation, multiple linear regression, and extended ARMA ( p , q ) models with exogenous climatic variables. The results showed that, according to the minimum of the mean absolute error, ARMA ( p , q ) models with optimally selected order parameters, and an adaptive window, had some advantages compared to other single-imputation methods. The ARMA ( p , q ) model produced a good quality of gap-ﬁlling in time series with the mean absolute error of 0.19 ◦ C and 0.08 Vol. % for soil temperature and moisture content, respectively. The ﬁndings supplemented the methodology of hydrological data processing and the development of digital tools for the online monitoring of climate and soil properties in agriculture.


Introduction
The monitoring of soil temperature and moisture at different depths is vital in agriculture to plan seasonal soil tillage operations, soil water conservation, and evaluation of crop growing conditions. Extreme soil temperature and low soil moisture negatively affect crop growth [1]. These soil properties significantly influence crop yield [2][3][4] and can cause physical damage and yield loss. Moreover, the findings within interdisciplinary studies [5] confirmed that the amount of moisture in the soil is linearly related to microbial respiration and affects their interaction. The data on soil temperature and soil moisture help in estimating the soil heat flux, which can provide data for the analysis of water balance and actual evapotranspiration (ETa) [6][7][8], with water content being the primary indicator to assess crop water stress [6,9,10]. Finally, soil temperature and moisture monitoring are in short-term gaps. However, other studies [28,29] prefer simple linear interpolation or single-imputation techniques, using mean or median values to fill short-term gaps (<5 h). What quality of gap-filling, based on such a convenient approach, do we receive? What is the share of lost information that could have been recovered by using more sophisticated methods? According to Y. Gao, C. Merz, G. Lischeid, and M. Schneider [18], families of adaptive statistical models-ARIMA (an autoregressive integrated moving average) model and ARCH (autoregressive conditional heteroscedasticity) model-are the most suitable for the nature of hydrological data. The application of these methods depends on the stationary of the original stochastic process. Based on this expert knowledge, the main research question was-What level of performance of the ARIMA models do we have in gap-filling in hourly time series of the soil temperature and soil moisture of topsoil, taking into account climatic settings? Despite a wealth of studies in soil hydrology, there is a shortage of studies that focused on the comparative analysis of methods performance to estimate missing values in hourly time series of soil temperature and soil moisture, measured by a weighable gravitational lysimeter station at topsoil.
The study objectivities are (1) to test extended ARIMA models with exogenous climate variables as statistical tools for gap filling in hourly time series of soil temperature and soil moisture at topsoil, (2) to compare the accuracy of missing values recovering based on three methods (linear interpolation, multiple linear regression and extended ARIMA model), and (3) to make recommendations regarding gap-filling techniques depending on the duration of a data gap.

The Monitoring Network
The network for monitoring climatic and soil hydrological settings was established in the dry steppe of the Kulunda Plain. This network was installed on the private farm "Partner" located in the village of Poluyamki in the Mikhaylovsky district in the Altai Krai. The installation was undertaken in the period 2012-2013 as a part of the Kulunda Project [30]. A weather station (WS) was installed, manufactured by German company "ecoTech" and a weighable gravitation lysimeter station (LYS) was also built, manufactured by "UGT-Muencheberg" (Germany). The equipment design and its parameters were already described in detail in several papers [8,12,28]. Two soil monoliths were extracted from a cultivated field (Monolith 1) and grassland (Monolith 2) separately to install LYS. The grassland had not been cultivated since 1950, while the arable land was cultivated according to the conventional tillage practice. Each monolith had a 1 m square of surface area and 2 m depth. According to FAO guidelines [28], the soils in monoliths belonged to Calcic Chernozems. The time of data registration was synchronized at both stations.
The equipment configuration allowed for sending data online through a regional mobile communication network and storing them in a data server at the digital platform of the Altai State University. The monitoring network continuously recorders the climate and soil properties in the studied area. This allows us to analyze and predict crop growing conditions in a short-term period. Additionally, we can detect the crop water stress and days of droughts, analyze water balance and actual evapotranspiration.
The equipment continuously measured climatic and soil hydrological indicators in a minute interval and automatically computed them to hourly data. The WS measured the following meteorological data-air temperature ( • C), air humidity (%), atmospheric pressure (hPa), solar radiation (W/m2), wind speed (m/sec), wind direction (decimal degree), and precipitation (mm). The lysimeter measured temperature ( • C), moisture (Vol. %), and water retention (hPa) at three different depths of soil (30,50, and 120 cm).

Available Data and Description
We focused on soil temperature and moisture content data, measured by the lysimeter station at 30 cm in the arable land (Monolith 1) and natural steppe vegetation (Monolith 2). The most critical period for recovering the lost data is the crop growing season, which Land 2021, 10, 579 4 of 17 usually starts in May and ends in September in the studied area. To analyze the performance of different gap-filling routines, a sample of artificial data gaps of various durations from short (1-12 h) and middle (13-24 h) to long (25-50 h) had to be investigated. It demanded a complete time series of meteorological and lysimeter measurements. We obtained a complete hourly time series of climate data for the whole period (2013-2019). The lysimeter data had some considerable gaps, which significantly limited the available data for our analysis. We used the entire time series from August to September in 2013, June 2014, from May to September in 2015 and 2016, and from May to August in 2017, 2018, and 2019. Thus, the test period covered different seasons of crop growing and climatic conditions.
In general, data preparation techniques for WS and LYS were similar. First, we converted hourly cumulative data of precipitation measured through multisensory "Vaisala" (WS) to absolute values per hour. Second, all data during system error or noticeable outliers were removed. Finally, we shaped a cross-sectional dataset consisting of the complete hourly time series of meteorological and lysimeter data (soil temperature and moisture content) for the vegetation periods of 2013-2019 (14,288 observations, in total for each monolith). The climate conditions and features of soil measurements at 30 cm are summarized in Appendix A.
According to the Shapiro-Wilk normality test [31], soil temperature and moisture content observations were not normally distributed (p-values < 0.001). The yearly histograms (Appendices B and C) also illustrated a significant asymmetry of probability distributions. Preliminary analysis of relationships between climate conditions and two studied soil properties in the monoliths 1 and 2 showed a moderate correlation between air temperature, atmospheric pressure, and the soil temperature at 30 cm (Table 1). The soil moisture content showed no correlation with all climatic indicators. Due to the problem of stationarity and heteroscedasticity, which was raised by Y. Gao, C. Merz, G. Lischeid, and M. Schneider [18], we tested the hypothesis of a stationary time series of soil properties based on the Augmented Dickey-Fuller Test (ADF) [32]. The test confirmed that both time series for each year were stationary (p-values for the null hypothesis "X is a non-stationary time series" were found less than 0.1). Therefore, considering the revealed statistical properties, we expected the extended ARIMA (p,0,q) models with exogenous climate variables to be suitable for the studied time series.
The following simulated the artificial data gaps-parts of the given time series were randomly dropped with a variation of a gap length from 1 to 50 h. We repeated data gap simulation for each dataset separately (Monolith 1 and Monolith 2) and obtained two samples of gap cases that had different duration varying from short (1-12 h) and middle (13-24 h) to long (25-50 h). In total, we tested 188 gaps, 93 gaps were generated in the time series of soil temperature and 95-in the time series of soil moisture. Simulating data gaps, we supposed that the interval between adjacent data gaps must be more than four days because previous observations have to be used to train the gap-filling routine. Thus, we obtained two random samples of gaps in the time series of soil temperature and moisture for monoliths 1 and 2.

Methods
The basic idea behind this analysis is to estimate missing values according to a method in each data gap using a dataset within a window covered K complete observations around a data gap. We also supposed that using climatic information can significantly improve gap-filling results. In the first stage, we trained models and estimated missing values for each data gap separately. In the second stage, the idea was to compare the estimated missing values with observed values using the original data and output the goodness of fit indicators.
In this regard, we introduced some notation. Let y(t) be an estimated value within a data gap and t-an index of the calculation period, ordered from 1 to T (the indexes of the first and the last missing values in a data gap). According to this notation, if the observed values y obs (t) are known, the bias error of treated missing values for each t in a data gap can be calculated as We tested the following different methods: (Method 0) simple linear interpolation, (Method 1) multiple linear regression of climate variables, and (Method 2) different variants of extended ARIMA(p,0,q) models with exogenous climate variables.
Method 0. The method of linear interpolation ("zero model"-M0) means simple filling of each data gap with constant growth per hour: where a-a parameter computed as a = (y(T + 1) − y(0))/(T + 1), y(0) and y(T + 1)-observed values at the beginning and the end a gap, respectively.
Replacing missing values according to Method 0 was considered as a simple and more convenient procedure, which, as expected, gave the worse quality of missing data recovering. However, the accuracy of this method was used to evaluate how other methods performed gap-filling better in the studied time series.
Regression-based imputation methods suppose the replacement of missing values by predicted values from regression estimation. The advantage of these methods is the exploitation of all observations, particularly meteorological data and, in the case of ARIMA models, previous observations to estimate values that are expected under given conditions. Method 1. The generalized linear model (GLM) recovers missing values based on observed climate variables C(t) in a given period (t = 1, 2, . . . , T) with a time-of-day component i(t): where C(t) = (C 1 (t),C 2 (t), . . . ,C 6 (t))-a vector of independent climate variables including air temperature ( • C), air humidity (%), atmospheric pressure (hPa), solar radiation (W/m 2 ), wind speed (m/s), and rainfall (mm); i(t)-a dummy variable associated with the timeof-day; a 0 , a 1 , . . . ,a 6 , β-parameters reflected the effect of independent variables on the outcome variable, computed according to maximum likelihood estimation (MLE). Due to the revealed stationarity of the time series of soil temperature and moisture, the ARMA model notation as a case of ARIMA (p,0,q) is more suitable to generate values to fill data gaps [18,30,33].
Method 2. The linear autoregressive moving average ARMA(p,q) model with exogenous variables is commonly specified as: where ρ(L p ) = 1 − ρ 1 L − ρ 2 L 2 − . . . − ρ p L p -the autoregression component AR(p) with an order parameter p, ϕ(L q ) = 1 + ϕ 1 L+ ϕ 2 L 2 + . . . + ϕ q L q -the moving average component MA(q) with an order parameter q, L p and L q are the lag operators with orders of p and q, respectively, in shorthand notation, where L j = y(t − j) indicates the value of y observed at j hours before y(t), a C(t)-linear regression component of independent climate variables notated according to the equation (3). We tested three variants of the ARMA model: AR(p), MA(q), ARMA(p,q) with adaptive parameters p and q. The orders p and q are given in units of time in a time series. The parameters p and q were optimally selected according to the minimum value of an information criterion Akaike (AIC) from a set of possible orders (p = 0, . . . ,p max ; q = 0, . . . ,q max ).
For each data gap in the studied time series, all models were trained separately. To estimate missing values according to the linear interpolation method (2), we used the start point and the endpoint of a data gap. To learn regression models (3)-(4), training datasets were extracted from the original time series, which had artificial gaps. Each training dataset represents parts of the time series within a window covered K observations before the beginning of a data gap and K observations after and consisted of all climatic variables and observed values of a dependent variable. To identify coefficients of GLM (3), we used the window of 2K observations, as shown in Figure 1. Taking into consideration the specifics of ARMA models, we trained them using K-hours from previous data. Estimation of missing values in a data gap was based on the identified models (3)-(4) using given climatic data for the window size K. We varied K from 2 to 4 days (48 to 96 h) for the models (3) and (4) and found the optimal window size according to the minimum root mean squared error. Figure 1 shows the main principles of trained dataset selection for all used methods.
where ρ(L p ) = 1 − ρ1 L − ρ2 L 2 − … − ρpL p -the autoregression component AR(p) with an order parameter p, φ(L q ) = 1 + φ1 L+ φ2 L 2 + … + φqL q -the moving average component MA(q) with an order parameter q, L p and L q are the lag operators with orders of p and q, respectively, in shorthand notation, where L j = y(t−j) indicates the value of y observed at j hours before y(t), a C(t)-linear regression component of independent climate variables notated according to the equation (3).
We tested three variants of the ARMA model: AR(p), MA(q), ARMA(p,q) with adaptive parameters p and q. The orders p and q are given in units of time in a time series. The parameters p and q were optimally selected according to the minimum value of an information criterion Akaike (AIC) from a set of possible orders (p = 0,…,pmax; q = 0,…,qmax).
For each data gap in the studied time series, all models were trained separately. To estimate missing values according to the linear interpolation method (2), we used the start point and the endpoint of a data gap. To learn regression models (3)-(4), training datasets were extracted from the original time series, which had artificial gaps. Each training dataset represents parts of the time series within a window covered K observations before the beginning of a data gap and K observations after and consisted of all climatic variables and observed values of a dependent variable. To identify coefficients of GLM (3), we used the window of 2K observations, as shown in Figure 1. Taking into consideration the specifics of ARMA models, we trained them using K-hours from previous data. Estimation of missing values in a data gap was based on the identified models (3)-(4) using given climatic data for the window size K. We varied K from 2 to 4 days (48 to 96 h) for the models (3) and (4) and found the optimal window size according to the minimum root mean squared error. Figure 1 shows the main principles of trained dataset selection for all used methods. To compare the performances of modeling procedures for each data gap, we also extracted control datasets from the original data. We used control datasets to estimate the goodness of fit statistics (5)-(7) for each gap in the time series separately.
The mean bias error (MBE): To compare the performances of modeling procedures for each data gap, we also extracted control datasets from the original data. We used control datasets to estimate the goodness of fit statistics (5)-(7) for each gap in the time series separately.
The mean bias error (MBE): The mean absolute error (MAE): Land 2021, 10, 579 7 of 17 The root mean squared error (RMSE): The mean absolute percentage error (MAPE): All calculations were performed using R with packages "forecast", "aTSa" and "lubridate". To estimate coefficients of ARIMA (p,0,q) models and select optimal order parameters p and q, we used the auto.arima() function from a package "forecast" with stepwise selection option (stepwise = TRUE), which allowed us to exclude insignificant independent variables.

Results and Discussion
Following the data gap-filling procedure, we estimated coefficients in regression models (3)-(4) and obtained predicted values for each missing value. In particular, the computation results revealed that the input of the time-of-day variables and climate information significantly improved regression models. The statistical significance of climate variables varied depending on a gap length and climate conditions. The hourly sum of precipitation and the wind speed were insignificant factors to predict soil temperature at topsoil. Thus, these variables were dropped out from the specification of this model. Despite no correlation between soil moisture and climate data, some of the climate variables were statistically significant in soil moisture regression models. The exception was the precipitation, which was excluded from all regression models. The list of significant variables was different for gaps and was not of interest to our study. As preliminary calculations showed, the maximum values of optimal order parameters p and q for AR (p), MA (q), and ARMA (p,q) models were 5 and 3, respectively. Therefore, to minimize computation time, optimization of order parameters was conducted by varying p from 0 to 5 and q from 0 to 3. In total, for 188 data gaps, we identified 188 GLM models and 564 ARMA models. For each model, we found the optimal window size (K); and for ARMA models, optimal values of p and/or q were computed, as well.

Mean Error Analysis
To compare gap-filling techniques, we calculated the means of the bias error indicators (MBE, MAE, RMSE, and MAPE) for all methods ( Table 2).
The results showed that ARMA (p,q) model had significant advantages in gap-filling in the time series of the soil temperature at 30 cm. The mean bias error produced by the extended ARMA model ranged from −0.70 to 1.02 • C and had a mean value of −0.02 • C with a standard deviation of 0.19. Compared to multiple linear regression and simple linear interpolation, the accuracy of missing values estimations based on the ARMA model was significantly better. For example, MBE for linear regression varied from −2.1 to 1.94 • C with a mean value of 0.04 • C and a standard deviation of 0.47 • C. The minimum of the mean absolute error (0.19 • C) was also achieved using the ARMA model, while MAEs generated by linear interpolation and linear regression were about two times higher. According to MAPE, the lost information in data recovery decreased from 2.4% in the linear interpolation to 1.3% on average for the ARMA model. However, the maximum values of MAPE reached 10.7%, 14.4%, 15.0%, 16.6%, and 10.1% for the linear interpolation, multiple linear regression, AR (p), MA (q), and ARMA (p,q) models, respectively. The The mode values for the optimal orders of ARMA (p,q) were 2 and 0 for p and q in both soil temperature and moisture content models. The average optimal window size to get accurate prediction was 71 h of previous observations. These findings confirmed the high suitability of the ARMA model to the nature of soil temperature and moisture content dynamics.

Error Analysis in Gap Size Class
To compare the performance of methods in relation to gap size classes, we estimated MAEs and their standard deviations for each gap length ( Table 3). The results showed the different power of the models in missing value estimation.
The best technique of gap filling in soil temperature time series, as mentioned above, was also the ARMA model with the lowest MAEs and standard deviations for all gap classes. MAE for short gaps (<13 h) was 0.06 • C with a standard deviation of 0.06 • C. The missing values in middle-length gaps can be estimated with an expected error of 0.16 • C with a standard deviation of 0.13 • C. Long-length gaps can be filled with an expected accuracy of 0.26 • C and its standard deviation 0.19 • C. These values were significantly smaller than the standard deviations of hourly soil temperature in given months. MAEs decreased by half in extended ARMA models due to adding climate information. We should also note that the adaption of order parameters in the ARMA model played a crucial role in improving gap-filling results. The filling of short-length gaps in soil moisture time series produced better results via simple AR (p) and ARMA (p,q) models, according to MAE-0.04 and 0.03 Vol.%, respectively. These models had moderate advantages in filling data gaps due to lower MAEs and standard deviations. Estimating missing values in middle-length gaps was done with an expected accuracy of 0.07 Vol. % and the standard deviation of 0.06-0.07 Vol. %. The input of climate information did not significantly decrease the error for long-length gaps (25-50 h), which reached 0.11 Vol. % on average. Figure 2 presents the general dependence of RMSEs on data gap lengths. Linear interpolation, GLM, MA (q), and AR (p) models produced higher RMSE with high variance compared to ARMA models. Importantly, ARMA models were generally better than other methods for all gap durations for soil temperature time series. Regarding soil moisture, we also observed an increasing RMSE, depending on the gap length and the advantage of ARMA models. should also note that the adaption of order parameters in the ARMA model played a crucial role in improving gap-filling results. The filling of short-length gaps in soil moisture time series produced better results via simple AR (p) and ARMA (p,q) models, according to MAE-0.04 and 0.03 Vol.%, respectively. These models had moderate advantages in filling data gaps due to lower MAEs and standard deviations. Estimating missing values in middle-length gaps was done with an expected accuracy of 0.07 Vol. % and the standard deviation of 0.06-0.07 Vol. %. The input of climate information did not significantly decrease the error for longlength gaps (25-50 h), which reached 0.11 Vol. % on average.  Figure 2 presents the general dependence of RMSEs on data gap lengths. Linear interpolation, GLM, MA (q), and AR (p) models produced higher RMSE with high variance compared to ARMA models. Importantly, ARMA models were generally better than other methods for all gap durations for soil temperature time series. Regarding soil moisture, we also observed an increasing RMSE, depending on the gap length and the advantage of ARMA models.
(a) Soil temperature (b) Soil moisture content  Comparisons of missing value estimations for two studied land use systems showed insignificant differences in the mean bias error between arable land (Monolith 1) and natural steppe vegetation (Monolith 2). p-values of t-test comparing two means MBE for the null hypothesis were more than 0.1 for all methods of gap-filling. Nonetheless, the analysis of MAPEs showed moderate differences in accuracy of estimations (Table 4).

Examples of Gap Filling and Selection of the Best Model
To illustrate the effectiveness of the gap-filling routines, two examples of missing value estimation are presented. Figure 3 presents the examples of missing value estimations for a long-length gap (48 h) in time series of the soil temperature. According to a common approach, the best model has to give estimations closest to actual values. In the present case, the best filling model was ARMA (4,1), which produced MAE equal to 0.06 • C and MAPE-0.31%. Comparisons of missing value estimations for two studied land use systems showed insignificant differences in the mean bias error between arable land (Monolith 1) and natural steppe vegetation (Monolith 2). p-values of t-test comparing two means MBE for the null hypothesis were more than 0.1 for all methods of gap-filling. Nonetheless, the analysis of MAPEs showed moderate differences in accuracy of estimations (Table 4).

Examples of Gap Filling and Selection of the Best Model
To illustrate the effectiveness of the gap-filling routines, two examples of missing value estimation are presented. Figure 3 presents the examples of missing value estimations for a long-length gap (48 h) in time series of the soil temperature. According to a common approach, the best model has to give estimations closest to actual values. In the present case, the best filling model was ARMA (4,1), which produced MAE equal to 0.06 °C and MAPE-0.31%.  As shown in Figure 4, the regression estimation of missing values in soil moisture time series had a high deviation, but the best models were ARMA (1,2) produced MAE-0.095 Vol. % and MAPE-0.31%. As shown in Figure 4, the regression estimation of missing values in soil moisture time series had a high deviation, but the best models were ARMA (1,2) produced MAE-0.095 Vol. % and MAPE-0.31%.    The results showed that MA (q) models were the least effective in gap-filling. The multiple linear regression models gave us low accuracy of estimations, even with the timeof-day variables, and can hardly be considered as a good tool of data gap filling in soil hydrological measurements. As was noted by Gao Y., Merz C., Lischeid G., and Schneider M. [18], regression methods tend to overestimate correlations and decrease variance. Moreover, according to a similar study [34], gap-filling in hourly time series of soil moisture showed the highest RMSE compared to other techniques. Simple techniques, such as replacement by monthly average, showed significantly better results. Linear interpolation cannot principally describe dynamics of soil hydrological settings during the day, which is important for using such estimations, especially if we face long-length gaps. Notably, the effectiveness of gap-filling methods with adaptive parameters such as ARMA (p,q) is limited if serial gaps become frequent and happen more often than 72 h [35].     The results showed that MA (q) models were the least effective in gap-filling. The multiple linear regression models gave us low accuracy of estimations, even with the timeof-day variables, and can hardly be considered as a good tool of data gap filling in soil hydrological measurements. As was noted by Gao Y., Merz C., Lischeid G., and Schneider M. [18], regression methods tend to overestimate correlations and decrease variance. Moreover, according to a similar study [34], gap-filling in hourly time series of soil moisture showed the highest RMSE compared to other techniques. Simple techniques, such as replacement by monthly average, showed significantly better results. Linear interpolation cannot principally describe dynamics of soil hydrological settings during the day, which is important for using such estimations, especially if we face long-length gaps. Notably, the effectiveness of gap-filling methods with adaptive parameters such as ARMA (p,q) is limited if serial gaps become frequent and happen more often than 72 h [35]. The results showed that MA (q) models were the least effective in gap-filling. The multiple linear regression models gave us low accuracy of estimations, even with the time-of-day variables, and can hardly be considered as a good tool of data gap filling in soil hydrological measurements. As was noted by Gao Y., Merz C., Lischeid G., and Schneider M. [18], regression methods tend to overestimate correlations and decrease variance. Moreover, according to a similar study [34], gap-filling in hourly time series of soil moisture showed the highest RMSE compared to other techniques. Simple techniques, such as replacement by monthly average, showed significantly better results. Linear interpolation cannot principally describe dynamics of soil hydrological settings during the day, which is important for using such estimations, especially if we face long-length gaps. Notably, the effectiveness of gap-filling methods with adaptive parameters such as ARMA (p,q) is limited if serial gaps become frequent and happen more often than 72 h [35].

Conclusions
Missing data comprise a common problem in meteorological and hydrological studies, which causes substantial difficulties for most statistical models, decreases the effectiveness of introducing online digital tools in agriculture. In this study, we compared various gapfilling methods in the time series of topsoil properties measured by the lysimeter station for two cases-arable land and natural steppe vegetation during vegetation periods of 2013-2019. Here, we applied different methods from conventional multiple linear regression practices and linear interpolation to the more sophisticated method-the extended ARMA model with exogenous climate variables. Despite the relatively low MAEs of gap-filling and the simplicity of the method, linear interpolation cannot describe non-linear hourly dynamics of soil temperature and moisture, especially if a gap length is longer than 6 h. Multiple linear regression and MA (q) models showed lower accuracy compared to the extended ARMA model. We confirmed that the ARMA model with exogenous climate variables has some advantages for missing values estimation in both soil time series under different climatic conditions. First of all, the considered approach has benefits, due to adaptive parameters of autoregressive and moving average components (p,q), and the adaptive window (K). Moreover, the input of climate information into the ARMA model demonstrated a significant improvement. Therefore, the ARMA (p,q) model can be used as a good tool for the prediction of soil properties and their dynamics, based on known climatic data and a part of the previous observations. This is important in case of shortterm equipment failure. However, the ARMA(p,q) model has some disadvantages: first, demand for sufficient dataset to train a model; second, the dependence of results on a criterion of selecting optimal parameters p and q. Here, we used AIC to select the optimal order parameters. Further studies will be challenged with at least three problems: the first criterion justification to select the best ARMA(p,q) model; the second analysis of gap-filling methods performance using a wider set of climatic events from droughts to more humid conditions, and a winter season; the third comparisons of the gap-filling effectiveness using various ARIMA, ARCH and GARCH models, as well as neural networks.
Author Contributions: These authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.