## Appendix A Kriging Interpolation  
The theoretical bases of geostatistical methods have been described in detail by [PERSON] [73] and [PERSON] [74], among others. Here, we just summarize some main aspects to consider for performing a 2-D spatial interpolation by kriging. In geostatistical methods, the variable of interest, Z, is considered as a stochastic spatial random process, i.e., a random function of space \(\mathbf{x}\), Z(\(\mathbf{x}\)), where \(\mathbf{x}\) is the vector containing a pair of coordinates, such as longitude and latitude. In this framework, an observed value of the variable of interest at a sampling location, \(z(\mathbf{x_{i}})\), represents a realization of the random variable Z(\(\mathbf{x_{i}}\)). The set of random variables \(\{Z(\mathbf{x_{1}}),Z(\mathbf{x_{2}}),\dots,Z(\mathbf{x_{m}})\}\) constitutes a random function for which only one realization, given by the set of observed values at the \(m\) sampling locations \(\{z(\mathbf{x_{1}}),z(\mathbf{x_{2}}),\dots,z(\mathbf{x_{m}})\}\), is known. The purpose of kriging is to predict data values at ungauged locations \(\mathbf{x_{k}}\) using information available at the \(m\) sampling locations \(\{\mathbf{x_{1}},\mathbf{x_{2}},\dots,\mathbf{x_{m}}\}\). This can be carried out by expressing Z(\(\mathbf{x_{k}}\)) as a linear combination of \(\{Z(\mathbf{x_{1}}),Z(\mathbf{x_{2}}),\dots,Z(\mathbf{x_{m}})\}\) as follows:  
\[\hat{Z}(\mathbf{x_{k}})=\sum_{i=1}^{m}w_{i}Z(\mathbf{x_{i}}) \tag{1}\]  
where the weights \(w_{i}\), for i = 1,... m, are obtained by solving the kriging system that provide the optimal (minimum variance) unbiased linear estimator [73]. In particular, ordinary kriging assumes the hypothesis of constant unknown mean, \(\mu\), over the search neighborhood of \(\mathbf{x_{k}}\). The ordinary kriging system can be written as follows:  
\[\left\{ \begin{array}{l}\sum\limits_{j=1}^{m}w_{j}\gamma\big{(}\big{\|} \mathbf{x_{i}}-\mathbf{x_{j}}\big{\|}\big{)}-\mu=\gamma(\big{\|}\mathbf{x_{i}} -\mathbf{x}\big{\|})\text{ for i = 1, \dots, m}\\ \sum\limits_{j=1}^{m}w_{j}=1\end{array} \right. \tag{2}\]  
where \(\gamma\) is the semi-variogram of the random process/field. Under the hypothesis of an intrinsically stationary, isotropic random field Z(\(\mathbf{x}\)), the semi-variogram of the field unambiguously describes the spatial autocorrelation of the sampling locations and it can be estimated from data as follows:  
\[\hat{\gamma}(\mathrm{h})=\frac{1}{2}\mathrm{E}[Z(\mathbf{x}+h)-Z(\mathbf{x})]^ {2} \tag{3}\]  
where \(h\) is the varying distance (lag) between two sampling locations. Then, the experimental semi-variogram has to be fitted by a suitable parametric analytic function to define the theoretical semi-variogram of the field, \(\gamma\) that is needed in the system of Equation (2). Among various analytical function used to fit the experimental semi-variograms, here we recall the exponential curve:  
\[\gamma(h)=n+s\big{[}1-\exp\!\left(-\frac{h}{r}\right)\big{]} \tag{4}\]where \(h\) is the lag distance; \(n\), \(s\) and \(r\) are parameters of the curve called nugget, partial sill and range, respectively. The nugget is the value of the semi-variogram at lag equal to zero and it is theoretically zero. However, in presence of small-scale variability and/or in presence of measurements errors, the nugget is not zero (nugget effect). The range defines the distance within sample locations are spatially correlated while the value that the semi-variogram model reaches at a distance equal to the range is called the sill. The partial sill is the sill minus the nugget. Appendix B The Network of Temperature Sensors and Sensitivity of the Kriging Estimator to Network Density  
Figure 11 shows, on the map of the region, the location of the 59 ground stations that are equipped with air temperature sensors along with the 18 ground automatic weather stations (AWSs) equipped with a complete set of weather sensors.