Improvement of Spatial Interpolation of Precipitation Distribution Using Cokriging Incorporating Rain-Gauge and Satellite (SMOS) Soil Moisture Data

Precipitation data provide a crucial input for examining hydrological issues, including watershed management and mitigation of the effects of floods, drought, and landslides. However, they are collected frequently from the scarce and often insufficient network of ground-based rain-gauge stations to generate continuous precipitation maps. Recently, precipitation maps derived from satellite data have not been sufficiently linked to ground-based rain gauges and satellite-derived soil moisture to improve the assessment of precipitation distribution using spatial statistics. Kriging methods are used to enhance the estimation of the spatial distribution of precipitations. The aim of this study was to assess two geostatistical methods, ordinary kriging (OK) and ordinary cokriging (OCK), and one deterministic method (i.e., inverse distance weighting (IDW)) for improved spatial interpolation of quarterly and monthly precipitations in Poland and near-border areas of the neighbouring countries (~325,000 or 800,000 km2). Quarterly precipitation data collected during a 5-year period (2010–2014) from 113–116 rain-gauge stations located in the study area were used. Additionally, monthly precipitations in the years 2014–2017 from over 400 rain-gauge stations located in Poland were used. The spatiotemporal data on soil moisture (SM) from the Soil Moisture and Ocean Salinity (SMOS) global satellite (launched in 2009) were used as an auxiliary variable in addition to precipitation for the OCK method. The predictive performance of the spatial distribution of precipitations was the best for OCK for all quarters, as indicated by the coefficient of determination (R2 = 0.944–0.992), and was less efficient (R2 = 0.039–0.634) for the OK and IDW methods. As for monthly precipitation, the performance of OCK was considerably higher than that of IDW and OK, similarly as with quarterly precipitation. The performance of all interpolation methods was better for monthly than for quarterly precipitations. The study indicates that SMOS data can be a valuable source of auxiliary data in the cokriging and/or other multivariate methods for better estimation of the spatial distribution of precipitations in various regions of the world.

The spatial distribution of rainfall is frequently assessed using ground-based raingauge stations [15,16] since they provide the most precise and reliable measured data [17]. However, the network of the stations varies considerably and is often scarcely and not adequately distributed to obtain sufficiently dense point measurements for generating continuous high-quality rainfall maps in target areas [7,18,19]. Hence, to gain a suitable spatial distribution of precipitations based on point data, spatial interpolation methods should be used [18,20,21]. The recent literature review [22] may suggest a conclusion that, due to the particularly complex spatiotemporal variability and physical mechanism of rainfall, acquisition of rainfall data of high quality and resolution is still challenging.
Several spatial interpolation methods employing rain-gauge records have been developed (e.g., conventional trend surface analysis Thiessen polygons, deterministic methods such as inverse distance weighted (IDW) estimating values at unsampled points by the weighted average at surrounding points [23,24], and more complex geostatistical methods (e.g., kriging) that count the spatial dependence between rain gauges and the spatial arrangement around the forecasting place) [25][26][27][28]. The suitability of each interpolator depends on the environmental conditions, and hence, a comparative study aiming at the selection of the best method in the target area is recommended [4,15,29].
Precipitations frequently exhibit high spatial variability related to topographic and climate factors [30][31][32]. Therefore, to improve the estimation of the spatial distribution of rainfall, a multi-element survey using the geostatistical cokriging method can be appropriate [5,15,18,29]. The cokriging method allows the joint use of, for example, a sparsely distributed variable such as rain gauges and a densely sampled secondary variable that complements the former [26,27]. A recent study conducted by Adhikary et al. [15] revealed that elevation employed as a secondary variable in the cokriging method significantly enhanced the estimation of rainfall in Australian catchments with mountainous and complex terrain. In another recent study carried out by Pellicone et al. [29] in a mountainous region of Calabria (Italy), ordinary cokriging (OCK) and kriging with external drift (KED) interpolators using elevation and distance to the coastline as secondary variables improved the spatial rainfall distribution. However, in a relatively flat region in Belgium, the OCK and KED methods with elevation as a secondary variable did not enhance the estimation of spatial rainfall distribution [33]. In Swiss conditions, cokriging using temporal rain-gauge data as a secondary variable improved the estimations of the precipitations compared with radar estimates [34]. In another study [35], blending satellite and ground precipitation observations using a Bayesian kriging approach significantly improved the satellite-derived estimates. The above literature suggests that further developments that include new secondary variables are required to improve the prediction of the spatial distribution of rainfall.
A secondary variable that has not been used so far in the cokriging method is the soil moisture content available from satellite data of the Soil Moisture and Ocean Salinity (SMOS) mission [36,37]. Soil moisture (SM), in general, rises upon the incidence of rain events [38,39] and affects precipitation events on the other hand [40,41]. The SMOS mission launched in November 2009 [42] offers soil moisture observations for the topsoil (up to several centimetres' depth) [43,44] using the interferometric radiometry method (1.4 GHz) from the orbit. SM data can be obtained every 1-3 days [45]. Validation has shown that SMOS SM data are well correlated with in situ soil moisture (R > 0.7) [46], exhibit reasonable root square mean error (RMSE) (≤0.04 m 3 m −3 ) [43,[46][47][48][49], and thus show good global performance [43,45]. Recent studies have demonstrated that assimilation of remotely sensed soil moisture has the potential to improve the quality of the near-real-time SMOSbased rainfall product [50,51] and to determine soil water resources [52].
It has been shown that the soil moisture-precipitation (SMP) feedback can be positive [53,54] or negative [55] depending on the environmental conditions. Several studies have demonstrated [40,53,56] that soil moisture information promotes rainfall generation, and the strength of this feedback is enhanced during warm seasons. In a study conducted by Vivoni et al. [40], the greater areal range of the storm maximum occurring in the high terrain of the Nacimiento Mountains was attributed to interaction between initial soil moisture and orographically induced rainfall. A negative SMP feedback was observed under a strong stability barrier at the top of the terrestrial boundary layer, which requires more sensible heat to induce appropriate turbulent mixing [55].
The main objective of this study was to examine whether satellite SMOS soil moisture data used along with rain-gauge records in the cokriging method are suitable for enhancing the estimation of the spatial distribution of quarterly (in the years 2010-2014) and monthly (in the years 2014-2017) precipitation over the area of Poland and near-border areas of neighbouring countries. The results of the cokriging method are compared with those of two methods (i.e., ordinary kriging and inverse distance weighting), which do not include satellite soil moisture data. Adequate evaluation of the spatial distribution of precipitation is important in a target area where rainfed agriculture is mostly practised and precipitation is the key uncertainty affecting water availability.

Study Area and Data Used
The study area included the whole territory of Poland and near-border areas of neighbouring countries (i.e., the Czech Republic, Slovakia, Ukraine, Belarus, Lithuania, and Russia), covering approximately 800,000 km 2 . The study area is dominated by soils with texture of loamy sands and loams. In the case of the examined area of Poland, it was approximately 325,000 km 2 . Quarterly precipitation data for a 5-year period (2010-2014) from 113-116 rain gauges retrieved from Tutiempo Network, S.L., Copyright 2018, located in a variety of networks were used ( Figure 1A). Monthly precipitation data were obtained from over 400 ground meteorological stations of the Institute of Meteorology and Water Management-National Research Institute (IMGW) ( Figure 1B). These free-access data are in ASCII format [57]. The original 1-day precipitation values were summed to monthly blocks in order to be consistent with the time scale of the SM from SMOS. It was assumed that a given weather station is representative of the entire SMOS pixel in which it is located. The rain-gauge stations were located at different altitudes varying from 1 to 1988 m a.s.l. The average altitude was 311 m with standard deviation as high as 320 m and, consequently, a coefficient of variation (average divided by standard deviation and multiplied by 100%) of 103%. The distribution of the altitudes exhibited a positive asymmetry and high slenderness. Skewness and kurtosis were 2.658 and 8.886, respectively. For this study, the SMOS L2 v. 650 datasets provided by the European Space Agency were examined. Based on the SMOS mission data processing algorithm, for each record assigned to a single DGG (Discrete Global Grid) node number, L2 retrieval has been carried out, under the condition that such pixel is not masked by no measurement value (−999 is assigned as an indicator of no measurements). The procedure for masking is defined by the series of quality flags defined in the mission Algorithm Theoretical Basis Document (ATBD) [58]. Additionally, variety of flags is applied to define the scene for which retrieval is conducted (including complex urban areas). The largest complications originate from external sources, such as radiofrequency interference (RFI) [58]. RFI contaminates the original signal, leading to exaggerated values of brightness temperature and, in consequence, unphysical values of SM. Thus, such pixels are masked by −999 (no measurement) value, before the L2 processing algorithm is applied. The datasets were downloaded from ftp://smos-ds-02.eo. esa.int/SMOS/L2SM/MIR_SMUDP2/. Next, further products were built from these data. The soil moisture contents in the topsoil (less than~10 cm) were extracted from the SMOS satellite data for about 5000 (for Poland and neighbouring countries) and 2000 (for Poland) points on the Discrete Global Grid (15 km grid) using the Icosahedral Snyder Equal Area (ISEA) map 4H9 projection. SMOS soil moisture was averaged quarterly in 2010-2014 [59] and monthly in 2014-2017. For the purpose of geostatistical analysis, a regularization procedure was performed (i.e., it was assumed that each pixel is represented by one point located in its centre).
We validated the soil moisture data based on measured soil moisture at the stations in the Podlasie and Polesie regions using classical statistics, Bland-Altman plot, concordance correlation coefficient, and total deviation index. The validation results were satisfactory [49].
Basic statistics including mean, minimal, and maximal values for annual and monthly rainfall and soil moisture were calculated in 2010-2014 and 2014-2017 ( Figure 2). As can be seen in Figure 2A Figure 2B). The coefficients of variation were similar in all study years (26.6-33.3%). The minimum and maximum soil moisture values were close to zero and 0.5 m 3 m −3 . The skewness (0.20-0.57) and kurtosis (0. 46-1.79) values indicate that soil moisture distribution was positively skewed and slightly narrow. The mean soil moisture in the individual quarters varied from 0.12 to 0.2 m 3 m −3 . The variability of soil moisture was, in general, highest in quarter I (to 40%) and lowest in quarter II (around 23%). The yearly standard deviations ranged from 0.04 to 0.051 m 3 m −3 . The standard variation was approximately 6.3% in quarter I and lower (up to ca. twofold) in the other quarters. The coefficient of variation was greater than 33% in each year and approximately 40% in quarter I of 2012 and 2013, with similar and lower variability in quarters II and III and increased variability in quarter IV. Histograms of soil moisture distributions indicate that the asymmetries were negative in quarter III and mostly positive in the other seasons (see Figure 3 in [59]). To obtain the normal distribution required in geostatistics, the variables with the highest asymmetry were transformed using the square root.
The  Figure 2D). The coefficients of variation were generally similar in all study years (23.9%-48.0%) except January 2017 (97%). The minimum and maximum of monthly soil moisture values were respectively close to zero and 0.748 m 3 m −3 (in wetlands). The values of skewness (−0.425 to 2.008) indicate that soil moisture distribution was slightly negatively and positively skewed, and those of kurtosis (−0.07 to 5.8) that it had a normal or narrow shape. Both the skewness and kurtosis values indicate that the monthly soil moisture distributions were close to the normal distribution.

Semivariograms and Cross-Semivariograms
The analysis of the spatial dependence and distribution of precipitation (z 1 ) and soil moisture (z 2 ) was performed using geostatistical methods. The normality of precipitation was obtained after square root transformation. After that, the soil moisture distribution was close to normal and thereby met the required condition of a stationary or quasi-stationary process. The experimental semivariogram γ(h) and cross-semivariogram between precipitation (z 1 ) and soil moisture (z 2 ) − γ 12 (h) for the distance h ( • ) were calculated from the following equations [60]: where N(h) is the number of pairs of points with values of [z 1 ( , distant by h, and x i is the spatial coordinate. For semivariograms and cross-semivariograms determined empirically, the following three mathematical models were selected: − spherical model: − exponential model: where γ(h) is the semivariance for internal distance class h, h is the lag interval, C 0 is the nugget variance ≥0, C is the structural variance ≥C 0 , A 0 is the range parameter, and C 0 + C is the sill. In the case of the spherical isotropic model, the effective range A = A 0 .
In the case of the exponential isotropic model, the effective range A = 3A 0 , which is the distance at which the sill (C + C 0 ) is within 5% of the asymptote. In the case of the Gaussian model, the effective range A = 3 0.5 A 0 , which is the distance at which the sill (C + C 0 ) is within 5% of the asymptote. In the case of the anisotropic model, the effective range , where A 1 is the range parameter for the major axis (φ) and A 2 is the range parameter for the minor axis (φ + 90). In the case of the exponential anisotropic model, the range (or effective range) is 3A 1 for the major axis and 3A 2 for the minor axis, φ is the angle of maximum variation, and θ is the angle between pairs. In the case of the Gaussian anisotropic model, the range (or effective range) is 3 0.5 A 1 for the major axis and 3 0.5 A 2 for the minor axis, φ is the angle of maximum variation, and θ is the angle between pairs. To evaluate anisotropy, the azimuth direction A z with the lowest semivariance values defined by smaller in the major direction (lower average semivariance) and largest in the minor (90 • -offset) direction was used.
The obtained mathematical functions of semivariograms and cross-semivariograms were used for the spatial analysis of autocorrelation or for the visualization, through estimation, of the rainfall value under consideration in space with the kriging or cokriging methods. In places where no samples had been taken, the data were estimated using the IDW (inverse distance weighting), OK (ordinary kriging), and OCK (ordinary cokriging) methods. Experimental semivariograms were calculated based on the rainfall dataset and soil moisture for each quarter and month with consideration of the effect of distribution anisotropy of both quantities on the semivariogram parameters. The spherical, exponential, and Gaussian functions were fitted (selected) to the empirical semivariograms for each quarterly and monthly dataset of rainfalls and soil moisture. Then, models (functions) with the minimum residual sum of squares and the highest R 2 values were selected.

Interpolation Methods
A deterministic inverse distance weighting (IDW) and two geostatistical methods (i.e., ordinary kriging (OK) and ordinary cokriging (OCK)) were used to assess the spatial distribution of precipitations.

Inverse Distance Weighting
Inverse distance weighting (IDW) estimates values of precipitations at unsampled points by the weighted average of observed data at surrounding points [23,24]. Estimation of values in places where no samples had been taken was made with the inverse distance weighting method using the following equation [60]: where z j *(h) is the estimated precipitation value at desired location j, z i is the measured sample value at point i, h j is the distance between z j *(h) and z i , s is the smoothing factor, and p is the weighting power.

Ordinary Kriging
Estimation of precipitation values in places where no samples had been taken was conducted with the kriging method [60]. The method yields the best nonbiased estimation of block values. This method also allowed for obtaining the minimum variance in the process of estimation. The estimator of kriging is a linear equation expressed by the following formula [60]: where N is the number of measurements, z(x i ) is the value measured at point x i , z*(x o ) is the estimated value at the point of estimation x o , and λ i is the weights. The weights are determined from a system of equations after inclusion of the condition of estimator nonbias and its effectiveness: Ordinary cokriging is a specific method for the analysis of random fields. It consists in the determination of covariance and reciprocal covariance as well as the crosssemivariogram function for specific parameters (i.e., precipitation (z 1 ) and soil moisture (z 2 )). The main advantage of the method is the possibility of indirect reconstruction of the spatial variability of precipitation, the measurement of which is difficult and expensive, through field analysis of soil moisture, which is easier to determine with standard measuring equipment or available from satellite observations.
During the estimation of values at sites where no samples had been taken, x o, can be made with the help of the estimation method known as the cokriging approach. The mathematical basis for cokriging is the theorem on the linear relationship of the unknown estimator z 2 * (x o ) expressed by the following formula [60]: where λ 1i and λ 2j are weights associated with z 1 and z 2 . N 1 and N 2 are the numbers of neighbours of z 1 and z 2 included in the estimation at point x o . Cokriging weights are determined from a system of equations with the inclusion of the condition of estimator nonbias and its effectiveness: where µ 1 and µ 2 are Lagrangian factors, and C 11 , C 12 , C 21 , and C 22 represent covariance between the variables. The relationships between the semivariance γ(h) and the covariance C are expressed by the equation In our study, the cokriging approach was used to enhance the estimation of spatial precipitation distribution using sparse data from rain-gauge stations and more densely sampled SMOS SM (as auxiliary variable) complementing the former. Cross-validation analysis is used for evaluating effective parameters for IDW, kriging, and cokriging interpolations. In this analysis, each measured point in the area is individually removed, and its value is estimated based on neighbouring measurement points. Then the point is replaced, and the next point is removed and estimated, and so on. Finally, the estimations are compared with measured values in all points, and statistical parameters are determined [60]. Table 1 presents the statistics of quarterly precipitations in 2010-2014. The mean precipitations were the lowest in quarters I and IV (80 and 150 mm) and considerably greater in quarters II and III (160 and 350 mm). The higher precipitations in the latter pair exhibited greater dispersion, as shown by the standard deviations. The coefficient of variability (CV) of the precipitations varied from about 30% to 62%. According to the literature [61], the CV values indicated moderate variability. The distribution of precipitations had a long and slender shape with positive asymmetries (values of 0.847-2.312) and large values of kurtosis (0.915-7.626). To meet the required condition of normality (or close to normality) in geostatistical analysis, the data were subjected to root square transformation, which led to lower asymmetries (0.218-1.327) and kurtosis (0.362-3.179).

Statistics of Rain-Gauge Data
The monthly precipitations in the years 2014-2017 varied largely from 0 to 438 mm, particularly during the summer (100-438 mm) ( Figure 2). The distribution of precipitations had a slight asymmetry and slightly blurred shape with positive asymmetries (values of 0.0-1.5) and small values of kurtosis (−0.2 to 2.1) and was close to the normal distribution.

Correlation Analysis
As can be seen from Table 2 It seems an obvious statement that soil moisture should increase when it rains and that there is no obstacle in reaching the soil surface. A positive and significant correlation between rainfall and soil moisture confirms the above statement for both quarterly and monthly rainfall data. The observed negative significant correlation between rainfall and soil moisture indicates that there are some factors in the studied areas that disturb this direction of soil moisture increase with increasing amounts of rainfall. Other factors changing this direction may be the soil texture and vegetation cover. We should also take into account the type of rainfall that occurs in a given area, whether it is frontal or convection rainfall limited to small areas. The type of rainfall will affect how quickly the water infiltrates into the soil, what part of the rainfall is captured by the vegetation, and how much water will evaporate from the soil and plant cover. Soil moisture in areas with coarse-textured soils can be less with more rainfall, and the inverse can be true in areas with fine-textured soils. Which of these factors and the relationships dominate in the study area and how they interact translate in turn into a correlation, whether it is significant, negative, positive, or insignificant.

Semivariogram and Cross-Variogram Models
The semivariograms and cross-semivariograms of paired precipitations with the SMOS satellite soil moisture for each year were calculated ( Table 3). The experimental values were best fitted (R 2 > 0.8) mostly to spherical, exponential, and Gaussian models. For quarterly precipitations and soil moisture, the best model to describe the spatial relationship was the exponential model. The ranges (A) of spatial dependence in semivariograms for quarterly precipitations were greater (from 1.26 • for III 2010 to 6.47 • for II 2010) with predominant values below 3 than for quarterly soil moisture contents (from 1.00 • for I 2010 to 4.08 • for IV 2013) with predominating values below 2 • (Table 3). In the case of cross-semivariograms, the ranges of spatial dependence between rainfall and soil moisture appreciably increased, reaching a maximum of 8.98 • with predominant values below 5 • . Anisotropy (A z ) varied for precipitations from 50 • to 120 • with predominant values above 100 • , for SMC from 0 • to 172 • with predominance below 120 • , and for P_SM from 12 • to 122 • with predominance above 100 • . The ratios C 0 /(C 0 + C) were in most cases (50 out of 60) <0.25, indicating strong spatial dependence; in the other cases, it was moderate [62] (0.25-0.75). It is worth noting that the spatial dependences (nugget/sill) were more frequently stronger in OCK (0.00-0.179) (except one, 0.403) than OK (0.018-0.25).
In the case of semivariograms for monthly precipitations and soil moisture, the best models to describe the spatial relationship were the spherical and the exponential models. As for monthly precipitation in the years 2014-2017, the effective ranges (A) of spatial dependence in semivariograms varied from 2.28 • for January 2015 to 5.77 • for October 2015 with predominant values from 3.73 • to 4.96 • . The corresponding range for soil moisture contents was from 2.49 for April 2015 to 5.84 for April 2014 with predominant values >4 • . In the case of cross-semivariograms, the best models to describe the spatial relationship were the Gaussian models. The effective ranges of spatial dependencies varied from 2.4 • to 4.79 • with predominant values above 3.61. The cross-semivariograms of monthly precipitations and soil moisture were negatively spatially dependent (10 times) and positively spatially dependent (6 times). Anisotropy (A z ) varied for precipitations and soil moisture from 58 • to 130 • with predominant values from 58 • to 95 • and P_SM from 100 • to 120 • . The ratios C 0 /(C 0 + C) were in cases (16 out of 32) <0.25, indicating strong spatial dependence; in the other cases, it was moderate (0.25-0.75). It is worth noting that the spatial dependences (nugget/sill) were strong in OCK (0.00-0.068). Table 3. Fitted semivariogram models for quarterly rainfall (P) and soil moisture (SM) data used in ordinary kriging interpolation method and cross-semivariogram models between rainfall and soil moisture (P_SM)-1 • in the ordinary cokriging method corresponds to about 100 km.

Comparison of the Interpolation Methods and Cross-Validation
The precipitation estimated with the use of the IDW, OK, and OCK methods was validated with the measured precipitation. The linear regression equation-y = ax + b, SEstandard error, R 2 -coefficient of determination, and SE Pre.-standard error prediction were calculated to compare the accuracy of each interpolator (Table 4). In the case of quarterly precipitation in the years 2010-2014, the ranges of regression coefficients (a), SE, SE Pre. (predicted), and R 2 for IDW and OK were 0. 552-1.428, 0.082-0.927, 29.1-108.9,  0.039-0.599 and 0.610-1.517, 0.072-0.339, 28.7-110.1, 0.040-0.634, respectively (Table 4). These values indicate a slightly better accuracy of OK than IDW. In turn, OCK indicates that the directional coefficients (a) of regression equations are above 1, up to 1.44, with small SEs (<0.032) and SE Pre. (<23) and with high values of R 2 (>0.94). The R 2 values for IDW and OK, which are the lowest in quarter I in most years (in 3 of 5 years), correspond with the lowest precipitation, but are similar in all the quarters in the case of OCK ( Table 4). The appreciably lower SE and SE Pre. values along with the higher R 2 for OCK than IDW and OK clearly indicate better performance of the former.
Irrespective of the type of interpolation method used, the values of R 2 were larger, and those of SE and SE Pre. were lower for monthly than for quarterly precipitations (Table 4).  Figure 3 presents example spatial distributions of quarterly precipitations derived from the OCK approach, which gave the best estimation among the three tested methods. They include the years 2010 and 2012 among the 5 years analysed (2010-2014). The precipitations (in all years) were lower than 450 mm in quarters I and IV and below 900 mm in quarters II and III. The precipitations recorded in the mountainous measuring points were even two times greater than the estimated values ( Figure 1) representing averages of both the high precipitation in the mountains and the lower precipitation in the neighbouring flat areas. In general, the variation in precipitation in the quarters is represented (characterized) by 3-4 and 10-13 colours, respectively, in the flat and less hilly areas with lower diversification of altitude and in the more mountainous areas with greater diversification ( Figure 1C). Generally, this differentiation was mostly latitudinal and was reflected in the greater values of anisotropy (about 100 • ) (0 • corresponds to anisotropy along the meridian) ( Table 3). A targeted growth in anisotropy from quarters II to III from the northwest to the southeast was observed in all years, except 2013, when precipitations of snow were abundant in quarter II (especially in April). The increases in precipitation from quarters IV to I (in all years) were not so pronounced or targeted. Figure 3 also presents the spatial distribution of monthly precipitations derived from the OCK method. Example distributions include the years 2014 and 2016 among the 4 years analysed (2014-2017). Variation of the monthly precipitation in Poland was generally smaller in winter and spring and much greater in summer and autumn (Figures 2 and 3). As in the case of quarterly precipitations, the orientation of monthly precipitations was generally mostly latitudinal and was reflected in the greater values of anisotropy with predominant values above 90 • . An anisotropy shift was observed from the northwest to the southeast to the west and the east from January to early October. Return to the previous orientation occurred in January of the following year, and the reorientation of anisotropy was repeated in the following years.

458
Interpolation methods are increasingly being used and improved as tools to enhance 459 estimation of the spatial precipitation distribution, as the availability of accurate obser-460 vation data is limited. This study showed that the spatial prediction of quarterly precip-461 itation data in Poland and near-border areas was more accurate when the ordinary 462 cokriging (OCK) method was used with rain-gauge data (P) as the main variable and the 463 satellite-derived soil moisture surface (SMOS SM) as an auxiliary variable (coefficients of 464 determination R 2

Discussion
Interpolation methods are increasingly being used and improved as tools to enhance estimation of the spatial precipitation distribution, as the availability of accurate observation data is limited. This study showed that the spatial prediction of quarterly precipitation data in Poland and near-border areas was more accurate when the ordinary cokriging (OCK) method was used with rain-gauge data (P) as the main variable and the satellite-derived soil moisture surface (SMOS SM) as an auxiliary variable (coefficients of determination R 2 = 0.944-0.992 between the measured and estimated precipitations), compared with both ordinary kriging (OK) (R 2 = 0.040-0.634) and inverse distance weighting (IDW) (R 2 = 0.039-0.599) using only the rain-gauge data. The greater suitability of the OCK method was supported by the smaller values of nugget and standard error prediction (SE Pre.) and the larger range of influence of the cross-semivariogram model than that of the direct (single) semivariogram for the above variables. As for monthly precipitation in the years 2014-2017, the performance of OCK (R 2 = 0.810-0.995) was considerably higher than that of IDW (R 2 = 0.352-0.844) and OK (R 2 = 0.348-0.865), similarly as with quarterly precipitation. The performance of all interpolation methods was better for monthly than for quarterly precipitations as shown by larger R 2 and lower SE and SE Pre. values.
Comparison of the present results and those from previous experiments indicates that the suitability of a given auxiliary variable depends on the topographic conditions of the studied areas. For example, in mountainous catchments of Australia at 25 to 1903 m a.s.l. [20] and the Calabria region in Italy [29] at <500 to 2266 m a.s.l., inclusion of elevation data as an auxiliary variable in the OCK method improved the prediction of the spatial distribution of orographically induced precipitations. However, this was not the case for the less hilly areas at 35 to 693 m a.s.l. in Belgium [33], where the OK and IDW methods performed better. Significant improvement of forecasting the spatial distribution of precipitations using OCK vs. OK and IDW observed in our study area with variable topography, including mountains, plateaus, plains, and valleys, may result from the fact that the topsoil (surface)-measured SMOS SM is sensitive to both small and large precipitation events in various topographic conditions. The suitability of satellite-derived SM data for improvement of prediction of spatial precipitation distribution based on rain-gauge data can be enhanced by their global availability [50,63,64]. The present analysis indicates that satellite-based remote sensing spatiotemporal soil moisture data can be a valuable source of an auxiliary variable for the cokriging and/or other multivariate (kriging) methods for better estimation of precipitations in various regions of the world. The performance of multivariate methods in the estimation of precipitation can be highlighted by recent attempts indicating that the spatial resolution of SMOS SM data can be increased by disaggregation of soil moisture from large to small pixels [65]. Additionally, it is possible to assess plant available water at a 0-50 cm depth based on the SMOS topsoil moisture using the Soil Water Index (SWI) and the exponential filter method as a proxy for the soil moisture in the root zone [47,63,66]. Furthermore, a combination of the previously developed approaches using satellite-derived soil moisture to improve satellite-based precipitation prediction [50,67] with the approach proposed in this study using satellite-derived soil moisture to improve rain-gauge-based precipitation prediction may have the considerable potential for upgrading precipitation prediction in future studies. Improvement of precipitation predictability by the combination of ground rain-gauge data with satellite-derived soil moisture data can be particularly advantageous in transitional climate zones under the influence of, for example, dry and wet climates [68], where they may reduce the strength of the positive soil moisture-precipitation feedback or cause even reversed feedback [53,69]. This may result from, for example, dissimilar evapotranspiration [70] and/or surface albedo [71] as well as the related net radiation influencing the amount and intensity of precipitation. In the area of the present study, the relatively low strength of the soil moisture-precipitation feedback can be affected by the transitional climate zone between the oceanic climate dominating in the north and west and the continental climate in the south and east of the area [72]. Additionally, the predictability of precipitations in the study area can be influenced by the spatial variability of the soil cover including coarse-textured permeable soils [73] that exhibit lower water holding capacity and heat capacity [74][75][76] compared with fine-textured soils. This effect can be associated with different warming-up and evaporation rates of the soils. Furthermore, coarse-textured soils display greater sensitivity to different precipitation amounts than other soils [14]. In connection with this, potential water deficits in the study area result mostly from the lack of water in the proper place [77] observed especially in the central part of the study area with soil cover dominated by coarse-textured soil [73]. Therefore, the accurate prediction of spatial distribution of precipitation obtained with the OCK method is anticipated to be beneficial in spatial planning and managing hydrological issues and soil water resources at various time frames [32].

Summary and Conclusions
The spatial distribution of quarterly (years 2010-2014) and monthly precipitations (years 2014-2017) was estimated in Poland and near-border zones of neighbouring countries using inverse distance weighting interpolation (IDW) and ordinary kriging (OK) based on precipitation data from rain-gauge stations. The methods were compared for the first time with ordinary cokriging (OCK) incorporating rain-gauge data (main variable) and soil moisture data from SMOS satellite measurements (auxiliary variable). In comparison with the IDW and OK methods, OCK was identified as a much better interpolator of the spatial precipitation distribution, as shown by the predominantly larger coefficients of determination (0.944-0.992 vs. 0.039-0.634) for quarterly precipitations than those for monthly precipitations (0.810-0.995 vs. 0.348-0.865) and by geostatistical measures including the nugget, standard error prediction, and range values. Therefore, we recommend the OCK approach for the generation of the most accurate (continuous) map of precipitations. The performance of all interpolation methods was better for monthly than for quarterly precipitations. This study emphasizes that combining rain-gauge precipitation data and satellite-based SMOS soil moisture products has a positive impact on the performance of spatial prediction of precipitations. The SMOS data products cover areas of different elevations and precipitations and thus partially reflect the topography effects on soil moisture. Since they are available globally, satellite-derived SM data have the potential to be used as an auxiliary variable in multivariate methods to enhance the estimation of precipitations in different regions of the world. Funding: This research was conducted (partially funded) under the project "Water in soil-satellite monitoring and improving the retention using biochar," no. BIO-STRATEG3/345940/7/NCBR/2017, which was financed by the Polish National Centre for Research and Development in the framework of "Environment, agriculture, and forestry"-BIOSTRATEG strategic R&D programme.
Institutional Review Board Statement: Not applicable.