## 2 Methods
### Gaussian Process with a Tapering Covariance Function for Local Variability
#### 2.3.4 Selection of Range Parameter \(\gamma\)  
The choice of the range parameter \(\gamma\) in the tapering covariance function influences the accuracy of the field estimation and the computational complexity of our proposed model. Given that maximizing the likelihood function in Equation (10) is to solve an unbiased estimating equation for \(\mathbf{\theta}\), it is suggested that an estimator of the sampling variability of \(\hat{\mathbf{\theta}}\) should be estimated using the robust information criterion in [33]. To achieve a trade-off between field estimation accuracy and computational complexity, we adopt the sampling variability (SV) and the model computational time (CT) as a criterion to select an appropriate \(\gamma\).  
Let \(\mathbf{U}(\mathbf{s},t;\mathbf{\theta})\) be an unbiased estimating function for \(\mathbf{\theta}\). The robust information criterion for \(\mathbf{U}\) is  
\[\mathcal{E}(\mathbf{U})=E\bigl{[}\dot{\mathbf{U}}\bigr{]}r\mathrm{E}[\mathbf{ U}\mathbf{U}\
u]^{-1}E\bigl{[}\dot{\mathbf{U}}\bigr{]}, \tag{12}\]where \(\dot{\mathbf{U}}\) is the matrix of derivatives of \(\mathbf{U}\) with respect to \(\boldsymbol{\theta}\),  
\[E[\dot{\mathbf{U}}]_{i,j}=-\frac{1}{2}\mathbf{tr}\left\{\left[ \frac{\partial\Sigma}{\partial\theta_{i}}\circ\mathbf{A}\right][\Sigma\circ \mathbf{A}]^{-1}\left[\frac{\partial\Sigma}{\partial\theta_{j}}\circ\mathbf{A }\right][\Sigma\circ\mathbf{A}]^{-1}\right\}, \tag{13}\] \[E[\mathbf{U}\mathbf{U}]_{i,j}=\frac{1}{2}\mathbf{tr}\left\{ \left[\left(\Sigma\circ\mathbf{A}\right)^{-1}\left[\frac{\partial\Sigma}{ \partial\theta_{i}}\circ\mathbf{A}\right][\Sigma\circ\mathbf{A}]^{-1}\right) \circ\mathbf{A}\right]\times\Sigma\left[\left(\Sigma\circ\mathbf{A}\right)^{-1 }\left[\frac{\partial\Sigma}{\partial\theta_{j}}\circ\mathbf{A}\right][\Sigma \circ\mathbf{A}]^{-1}\right)\circ\mathbf{A}\right]\Sigma\right\}. \tag{14}\]  
The sampling variability is defined by the diagonal elements of \(\mathcal{E}(\mathbf{U})^{-1}\), namely, the trace of \(\mathcal{E}(\mathbf{U})^{-1}\).  
The procedure of selecting \(\gamma\) is as follows (Figure 4): First, we obtain a pilot estimate \(\hat{\boldsymbol{\theta}}_{pilot}\). For example, we calculate the MLE for a small and randomly chosen subset of sensing observations. Second, we calculate the sampling variability \(SV(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i})=\mathbf{tr}[\mathcal{E}( \mathbf{G}(\mathbf{s},t;\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i}))^{-1}]\) and model computational time \(CT(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i})\) for an available value of \(\gamma_{i}\). The initial value \(\gamma_{1}\) is the minimal range parameter value. Third, we repeat the second step and calculate the sampling variability \(SV(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i})\) and model computational time \(CT(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i})\) for \(\gamma_{i}\) with \(i=2,\ 3,\ \ldots,\ n_{\gamma}\), where \(n_{\gamma}\) is the number of \(\gamma\); Finally, we select \(\gamma^{*}\) to construct a balance between the sampling variability \(SV(\gamma)=\left\{SV(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i}),i=1,\ 2,\ \ldots,\ n_{ \gamma}\right\}\) and the model computational time \(CT(\gamma)=\left\{CT(\hat{\boldsymbol{\theta}}_{pilot},\gamma_{i}),i=1,\ 2,\ \ldots,\ n_{ \gamma}\right\}\), that is, smaller sampling variability and shorter model computational time.