Stochastic Model of Spatial Fields of the Average Daily Wind Chill Index

: The objective of this paper was to construct a numerical stochastic model of the spatial field of the average daily wind chill index on an irregular grid defined by the location of the weather stations. It is shown in the paper that the field in question was heterogeneous and non-Gaussian. A stochastic model based on the real data collected at the weather stations located in West Siberia and on the method of the inverse distribution function that sufficiently well reproduce different characteristics of the real field of the average daily wind chill index is proposed in this paper. I also discussed several questions related to the simulation of the field on a regular grid. In the future, my intention is to transform the model proposed to a model of the conditional spatio-temporal field defined on a regular grid that allows one to forecast the wind chill index.


Introduction
During the last decades there arose many discussions about how to determine the thermal comfort and how to grade the thermal stress. These efforts have resulted in various models attempting to describe the thermal comfort and the resultant thermal stress (one of the first models was proposed by P.A. Siple and C.F. Passel in 1945 [1]). A large number of the so-called biometeorological indices that describe the combined thermal effect of various meteorological parameters (the wind speed, the air temperature and the relative humidity, the atmospheric pressure, etc.) on human beings have been proposed. In [2][3][4][5] one may find a comprehensive review of such biometeorological indices.
The most common approach to studying the properties of the time series and spatial/spatiotemporal fields of bioclimatic indices is the statistical approach (see, for example, [6][7][8][9][10]). The statistical approach could be easily transformed to a stochastic one. The obtained stochastic models may be used both to study the properties of a bioclimatic index and to forecast it. For example, in [11,12] two stochastic models of the bioclimatic index of severity of climatic regime are proposed and in [13][14][15] the same approach was used for the time series of the heat index.
In this paper, a stochastic approach was for the first time applied to studying the spatial fields of the average daily wind chill index (WCI). Earlier, this approach was successfully applied to simulation and forecasting of the WCI time series [16,17]. The field of the WCI was simulated on an irregular grid determined by the location of the weather stations. Several questions related to the simulation of the WCI spatial field on a regular grid are discussed. Both to determine the model input parameters and to verify the model proposed, the real data collected at the weather stations located in West Siberia were used.
There is a wide range of methods for simulation of Gaussian or non-Gaussian random processes and fields [18][19][20][21][22][23][24][25][26]. The problem is that most of the methods could be used only for simulation of some

The Wind Chill Index
The idea of the wind chill index was introduced by P.A. Siple and C.F. Passel in 1945 [1] after the United States Antarctic Expedition. Both the impromptu experiments that were conducted during the expedition (the results of these experiments have formed the basis of the WCI concept) and the proposed formulas for calculation of the WCI were criticized by many of the researchers (see, for example, [27,28]). Nevertheless, according to the [29], the WCI has become a standard meteorological term that is used to represent the hypothetical air temperature that would, under conditions of no wind, effect the same heat loss from unclothed human skin as does the actual combination of air temperature and wind speed. In [30], one may find some research background related to innovative formulas to estimate WCI. Different formulas for calculating the WCI and the restrictions under which these formulas could be applied are presented there.
In this paper, I used the definition of the WCI presented in [2,31,32] if the air temperature is below o 10 C and the wind speed is above 5km / h . In the case of the very low wind speeds, that is, less than 5km / h , the WCI is either not calculated at all or is defined according to the formula [33] 1 59 0 1345 5 If the air temperature is above o 10 C , the WCI is not calculated. In [2], one may find a table that presents the dependence of the health risk on the values of the WCI. The lower the WCI, the higher is the risk of frostbite and hypothermia. For example, if 27 10 W     , there is a low risk of frostbite and the hypothermia is possible if a person is outside for long periods without adequate protection. Meanwhile, if 47 40 W     , there is a high risk of frostbite: Exposed skin can freeze in 5-10 min.

The Real Data Description and Analysis
Both to determine the model parameters and to verify the model, the real data collected at the weather stations located in the southern part of West Siberia (Russia) from 1988 to 2017 were used [34,35]. I studied the properties of the average daily WCI and constructed the proper stochastic models for the two areas.  Figure 1 shows the positions of the weather stations in the areas in question.
It should be noted that the WCI is not measured at a weather station, but it could be easily calculated with the above-given definition, Equations (1) and (2), using the values of the air temperature and the wind obtained at the station. At all the weather stations considered, the air temperature and the wind speed data were collected eight times per day at the synoptic hours (00:00, 03:00, 06:00, 09:00, 12:00, 15  In order to construct a proper stochastic model, one should pay attention to the fact that the real fields W  of the WCI are heterogeneous. Figure 2 shows examples of the sample mean of the average daily WCI at the weather stations in the Area 1. One can easily see that the sample mean varies from one station to another. The properties of the spatial field of the average daily WCI are time-dependent also. This should be taken into account when one constructs a model of a spatio-temporal field. Since in this paper only spatial fields were studied, this time dependence is out of range of my interests.
The spatial fields of the average daily WCI were strongly correlated. Figure  respectively. The correlation function of the wind speed field decreased faster than the correlation of the WCI field, and the decrease rates of the correlation functions of the temperature and WCI fields were similar. To reduce the statistical error of the estimation of all the characteristics considered over a small sample of the real data, the sample size was increased artificially. To do it, the moving averaging procedure with a symmetric two-side smoothing window with a width of 2L 1  days was applied. To estimate the characteristic for the given day, data collected in that day, data collected during L days before the given day, and L days after the day in question were used. In this paper, L 2  .

Stochastic Model
The method of the inverse distribution function (MIDF) is used for the simulation of the random field of the average daily WCI [18,26,36]. This method allows one to simulate a non-Gaussian process with the given one-dimensional distributions and the given correlation structure. It is shown in [16] that the mixtures of the two Gaussian distributions will sufficiently approximate the sample histograms of the real average daily WCI. The correlation structure of the field W  is defined that is estimated on the real data. Let us recall the idea of the MIDF. In the framework of the MIDF, simulating of the field W  with the given densities  , 1, k g x k NS  and the given correlation matrix R comes to an algorithm (let us call it the Algorithm 1) with the three steps [11,16]: Step 1. At this step, a correlation matrix where the function is the distribution density of a bivariate Gaussian vector with zero mean, variance equal to 1 , the correlation coefficient  


Step 2. After calculating the matrix R' , a trajectory of the auxiliary standard Gaussian field

W '
 with the correlation matrix R' is simulated. This could be done using the Cholesky method or the spectral decomposition of the matrix R' (see [18]).


Step 3. Finally, the trajectory of the Gaussian vector W '  is transformed to the trajectory of the non-Gaussian vector W One should note that there is no analytical solution of the Equation (3). This is why one have to solve this equation numerically. Since the integral in the right part of Equation (3) is a continuous monotonically increasing function on [−1;1] (see, for example, [18]), the simplest method for solving the Equation (3) is the numerical inversion of the integral as a function of the correlation coefficient using its tabulated values. In this paper, the bisection method was used. If the matrix R , obtained at the first step using the Equation (3), is not positive definite, it must be regularized. In this paper, a method of regularization based on the substitution of negative eigenvalues of the matrix R with small positive numbers was used [18]. After the regularization, the matrix R required normalization. Steps 2 and 3 were repeated as many times as trajectories were required. For the transformation Equation (4) the tabulated values of Let me list the main advantages of the model proposed.


Many of the dynamic models that are used to simulate the fields of the air temperature and other meteorological parameters (using these models one may calculate the WCI) require as the input information some data related to the geographical properties of the area (whether the area is mountainous or plain, are there any huge water bodies or not, the elevation of the weather stations above the sea level, etc.). The stochastic approach does not require this additional information.  The most time-consuming step of the above-described algorithm is Step 1.
Step 2 and Step 3 do not require a lot of time to be conducted. Once the matrix R' is calculated, one may repeat Steps 2 and 3 as many times as needed to obtain the required accuracy of the estimations, and the increase in accuracy insignificantly influences the total computational time.


The model proposed could be easily transformed to a model of a conditional random field of the average daily WCI. For the time-series of the WCI, such transformation is presented in [17]. The conditional model allows one to forecast the bioclimatic index in question.
If one simulates the field of the average daily WCI only at the weather stations, the abovedescribed approach could be used (let us call it the Model 1). In the next section, the results of the verification of the Model 1 are given. However, if it is necessary to simulate the field not only at the weather stations but also in the nodes of some regular or irregular grid, one must determine the distribution density of the WCI in each node and determine the     are often used. Here Before one starts to simulate the field of the WCI on a grid with the correlation functions Equations (5) or (9), it is necessary to study how such an approximation influences the quality of the simulation of the field at the weather stations. Let the Model 2 and Model 3 be the models of the field of the WCI based on the MIDF, in which the matrix R was substituted with the matrices calculated with the Equations (5) and (9) with the parameters that minimize the functional Equation (10). In the next section, the results of the verification of the Model 2 and the Model 3 are given.

Verification of the Model
To verify a model, it was necessary to compare estimations of various characteristics that were based on the simulated and real data. Only such characteristics must be considered that, on the one hand, were reliably estimated by means of real data and, on the other hand, were not input parameters of the model. In this section, several examples of such characteristics used for verification of the Models 1-3 are given.
It should be noted that it is possible to simulate as many trajectories as needed to provide a required accuracy of the estimation when the simulated data are used. In this paper, for all estimations based on the simulated data, I attained the accuracy above 0.0001. Thus, in all the tables presented, the estimations based on the simulated trajectories are given with significant digits only.
From now on,  is a standard deviation of the characteristic under consideration when estimating with the real data. Since for the Models 1-3 only one-dimensional distributions of the field and its correlation structure were defined, it was not possible to write down theoretical formulas for the variance (and, therefore, for  ) of the characteristics of the field related to its multi-dimensional distributions. This is why the values  , presented in tables below, were numerically estimated.   Tables 3 and 4 show the examples of the corresponding estimations, obtained with the real and simulated data. When I considered the relatively high levels l , where the estimating of   p l was reliable (lines 1-5 in Tables 3 and 4 never happened. One could see that the Models 1-3 gave realistic results. This results were close to each other because the only difference between these models was the type of the correlation function (that characterizes the two-dimensional distributions of the field), and   p l was related to the NS  dimensional distribution of the WCI field. Another characteristic used for the comparison of the real and the simulated fields of the average daily WCI was the probability   , pr k l that at least at k stations the average daily WCI did not exceed the level l . Figure 5   The possible explanation is that the characteristic in question was closely related to the correlation structure of the field, and the correlation function Equation (9) used in the Model 3 was closer (in the sense of the functional Equation (10)) to the sample correlation coefficients of the real field of the average daily WCI. Figure 6 presents the estimations of   1 2 , , pd station station  for the pairs of stations Bakchar and Tomsk and Pervomajskoe and Bolotnoe. One can see that, despite the fact that all four stations are situated in the same climatic zone and the distance between Bakchar and Tomsk is almost equal to the distance between Pervomajskoe and Bolotnoe, there is a significant distinction between the decreased rates of   , , pd    as a function of  . The most probable explanation of this fact is that the correlation coefficient between the average daily WCI at the weather stations Pervomajskoe and Bolotnoe was higher than the correlation between the average daily WCI in Bakchar and Tomsk, and, therefore, the probability of the sufficient difference in values of the WCI at the first pair of stations was less than the corresponding probability at the second pair of stations.  AN l ), a fine grid with a small step in this argument was used for the verification. The results of this detailed numerical analysis showed that the trajectories obtained with the models proposed were close in their statistical properties to the real fields of the average daily wind chill index. Therefore, these models could be used for studying of those properties of the field that are unreliably estimated by means of real data. Since the Model 3 reproduced the characteristics of the real field a little more accurately than the Model 2, the Model 3 will be used for simulating of the spatial and spatio-temporal fields of the average daily wind chill index on a regular grid.

Conclusions
For the first time ever, the stochastic approach was used for studying and simulation of the random fields of the average daily wind chill index. It was shown in this paper that the spatial fields in question, defined as two areas located in West Siberia, were heterogeneous and non-Gaussian. I proposed a stochastic model that fairly well reproduces the statistical properties of the real spatial fields. This model could be used for estimating of those characteristics of the random field that could not be estimated from the real data in view of the small sample size of the real data.
In the future, our intention is to transform the proposed model to a model of the conditional spatio-temporal field defined on a regular grid that allows one to forecast the wind chill index. To achieve this goal, it is necessary to solve the next few problems. The first one is to determine such interpolation technique that allows one to interpolate values of the WCI from the weather stations to the grid nodes with the smallest error of the interpolation. The interpolation should let one take into account both spatial and temporal variability of the WCI. The next problem we have to solve is simulation of the nonstationary heterogeneous non-Gaussian field on a dense grid. Simulation of such a process on a grid with a large number of nodes is a time-consuming procedure and it requires also a lot of computer memory. It is quite probable that it will be necessary to substitute the Algorithm 1 with some approximate algorithm with the lesser computational complexity (see, for example, [23,24]). Finally, it would be necessary to transform the non-conditional model to a conditional one. This transformation is rather simple if the conditions are point conditions, when several values (either consecutive or inconsecutive) of the WCI are given. Since there is no efficient algorithm to simulate the conditional Gaussian time series with the interval conditions, the transformation of the nonconditional model to a conditional one is much trickier and requires further research.