Reconciling high resolution climate datasets using KrigR

There is an increasing need for high spatial and temporal resolution climate data for the wide community of researchers interested in climate change and its consequences. Currently, there is a large mismatch between the spatial resolutions of global climate model and reanalysis datasets (at best around 0.25o and 0.1o respectively) and the resolutions needed by many end-users of these datasets, which are typically on the scale of 30 arcseconds (~900m). This need for improved spatial resolution in climate datasets has motivated several groups to statistically downscale various combinations of observational or reanalysis datasets. However, the variety of downscaling methods and inputs used makes it difficult to reconcile the resultant differences between these high-resolution datasets. Here we make use of the KrigR R-package to statistically downscale the world-leading ERA5(-Land) reanalysis data using kriging. We show that kriging can accurately recover spatial heterogeneity of climate data given strong relationships with co-variates; that by preserving the uncertainty associated with the statistical downscaling, one can investigate and account for confidence in high-resolution climate data; and that the statistical uncertainty provided by KrigR can explain much of the difference between widely used high resolution climate datasets (CHELSA, TerraClimate, and WorldClim2) depending on variable, timescale, and region. This demonstrates the advantages of using KrigR to generate customized high spatial and/or temporal resolution climate data.


INTRODUCTION 25
Ongoing climate change is having wide-reaching effects worldwide. This has fuelled a need 26 for high spatial and temporal resolution climate data to be able to quantify and predict the 27 effects of climate change (Hewitt et  spatial and temporal resolution. However, due to the diversity of data sources and methods 35 used, all of these high-resolution datasets contain numerically different climate data, 36 particularly in locations with a low density of in-situ observations such as Alaska. This presents 37 a serious challenge for users of these datasets who need to know which dataset they can trust 38 for their particular purpose. Reconciling the differences between these different products is a 39 difficult task which is exacerbated by the lack of uncertainty metrics associated with the high-40 resolution data. None of the existing high-resolution climate products account for the 41 uncertainty in the underlying climate data, or in the downscaling technique, which can lead to 42 over-confidence of end-users in the validity of these products. One tool that has emerged which 43 may be used to address this challenge is KrigR (Kusch, Davy, in prep.). 44 45 KrigR is an R-Package which offers functionality for the retrieval and pre-processing of 46 ERA5(-Land) data as well as statistical downscaling of spatial products to high spatial 47

STATISTICAL INTERPOLATION 124
The KrigR package carries out statistical downscaling through kriging -a statistical 125 interpolation technique. Kriging is a two-step process that requires training data that we wish 126 to downscale, and co-variate data both at the resolution of the training data and at our target 127 spatial resolution (Chilès and Delfiner 2012). In the first step, we fit variograms to our training 128 data and establish covariance functions with our co-variate data at the training resolution. This 129 gives us functions which describe the spatial autocorrelation of our training data, and its 130 relationship with our chosen co-variate(s). During the second step we predict the value of our 131 variable at new locations, in this case at a higher spatial resolution, using co-variate data at the 132 target resolution. One major advantage to kriging is that it preserves the uncertainty obtained 133 when fitting the variogram, which gives us an uncertainty associated with the downscaled data. 134 In KrigR this uncertainty is given as a standard deviation of the uncertainty in the estimate 135 (Hiemstra et al. 2009). 136

137
There are a few important limitations to all statistical downscaling methodologies, and for some 138 variables there is simply no reasonable way to statistically downscale them. What cannot be 139 accounted for when using any statistical downscaling approach is the effect of dynamical 140 processes which occur only at the unresolved scales of the input data. For example, suppose 141 one wishes to downscale temperature data from a resolution of 100 km to a resolution of 1 km. 142 The dynamical model that was used to create the 100 km resolution data will account for how 143 air temperature varies with altitude, because this can be partially resolved at these scales, but 144 it will not include the effect of atmospheric circulation within a valley on air temperature. This 145 can become important when the atmosphere is stably stratified and cold pools (Mutiibwa et al. 146 2015) or frost pockets can form in topographical depressions. For this reason, dynamical 147 models of the atmosphere include descriptors of the un-resolved topography of the surface (e.g. heterogeneity of slope angle) which are used to account for such unresolved processes. 149 Whether or not statistical downscaling can account for these processes depends upon whether 150 the underlying process is represented in the training data, and whether the relevant co-variates 151 are used in the kriging. This essentially puts a limit on how large a change in resolution can be 152 accomplished using statistical downscaling. As a general guide we do not advise downscaling 153 to a resolution more than ten times finer than the training dataand so we have added a warning 154 in KrigR to the user to caution against this. 155 For some variables, such as precipitation, the processes that determine their spatial pattern at 156 finer resolutions than the training data are largely determined by atmospheric dynamics. 157 Therefore, no combination of topographical co-variates is going to enable us to statistically 158 downscale precipitation with high accuracy. We therefore do not recommend statistically 159 downscaling precipitation data. However, there can be alternatives which also tell us about the 160 water availability at high resolution, such as soil moisture, that we can successfully statistically 161 downscale by using the soil properties and topographical properties as co-variates as we have 162 done here. 163

DYNAMICAL AND STATISTICAL UNCERTAINTY 164
KrigR provides statistical uncertainty alongside downscaled products. This statistical 165 uncertainty is the uncertainty resulting from statistical downscaling and is given as a standard 166 deviation of the uncertainty in the estimate, σ Krig , at each timestep. The magnitude of this 167 uncertainty will depend upon (1) the robustness of the relationships between the target variable 168 and the co-variates, (2) the spatial variability of the training data, and (3) the change in 169 resolution between training and target resolutions (Chilès and Delfiner 2012). 170 In addition to this statistical uncertainty, the ERA5 reanalysis provides uncertainty due to the 171 dynamics of the climate system and the limited observations used to constrain the reanalysis, 172 henceforth referred to as dynamical uncertainty (Hersbach et al. 2020). In ERA5, this uncertainty information is provided as the 10 individual members of the ensemble used to 174 generate the reanalysis, and as a measure of the ensemble spread (standard deviation of the 10 175 members), which can both be acquired using KrigR. Here, we refer to this measure as is σ . Finally, to calculate the spread in the difference between the datasets at each timestep we first 195 calculated the differences between the datasets at each timestep e.g. for CHELSA this is given by: T' =TKrig-TCHELSA. We then fitted a normal distribution to these differences which gave us 197 both the mean difference between the datasets Err and the spread in the differences, σ , for 198 each gridcell. 199 to the surface temperature field as air temperatures tend to decrease strongly with altitude 206 (typically by 6.5 K km -1 ). Furthermore, the direction which sloping terrain is facing and its 207 steepness can also strongly affect surface air temperature by altering the amount of 208 downwelling solar radiation per unit area that is absorbed at the surface. Slope direction and 209 steepness can also strongly affect runoff, and hence soil moisture content. 210 Figure 1 shows the effect of using different co-variates when kriging surface air temperature 211 and soil moisture over the UK at a randomly chosen monthly time-step. We chose the UK as 212 an example due to large variations in topography as well as soil properties. Our downscaled 213 product shows that elevation is extremely important in determining the surface air temperature 214 and that the spatial pattern thereof closely matches that of elevation ( Figure 1A Including the slope steepness to the elevation-driven kriging ( Figure 1C) also changes the 220 estimate, but to a lesser degree with local changes of up to 0.1K. Doing so reveals clear patterns 221 in how the estimate changes, especially in the complex terrain of the Scottish Highlands: there 222 is a general warming at high altitudes and cooling in valleys. This pattern is due to the co-223 variance between elevation and slope steepness resulting in some component of the vertical 224 temperature gradient being assigned to steepness instead of elevation. By including slope 225 aspect with elevation ( Figure 1D), the estimate of surface air temperature changes by around 226 0.4K for select areas. However, in this case, there are less clear spatial patterns to how the 227 estimate changes as slope orientations can change completely from one grid-cell to the next. 228

RESULTS
For soil moisture, also shown in Figure 1 (E-H), the picture is very different. Accounting for 229 the soil properties ( Figure 1E) is important, but the addition of elevation ( Figure 1F) has little 230 effect on the downscaling output. Introducing slope steepness ( Figure 1G) also had little effect. 231 However, including slope aspect ( Figure 1H) as a co-variate can change the estimate by up to 232 0.1 kg kg -1 , a change of more than 15% from the estimate using soil properties alone. This tells 233 us that, on monthly-averaged timescales, soil moisture is strongly determined by the soil 234 thermal and hydrological properties. It is also more determined by the direction the terrain is 235 facing, and thus the amount of solar radiation it receives than by how steep the slope is or how 236 high the elevation. Note that the strength of these relationships is likely to vary with temporal 237 resolution, time of year, and geographical region. 238

DYNAMICAL UNCERTAINTY IS COMPARABLE TO STATISTICAL UNCERTAINTY 248
Uncertainties change with timescale. The magnitude of dynamical uncertainty can change from 249 hour to hour, but at hourly resolutions the dynamical uncertainty of surface air temperature can 250 be of the same magnitude or larger than the statistical uncertainty. Figure 2  can or should be used as a proxy for the other. 275

KRIGING RECOVERS SPATIAL VARIABILITY 276
Given the absence of independent, high resolution climate data with which to measure 277 downscaling skill of KrigR, we assess how well kriging can recover spatial variability from 278 coarse-grained ERA5-Land climate variables. First, we aggregate the ERA5-Land variable to 279 a coarser resolution (0.4 o ; 40km), then we use KrigR to statistically downscale the data back 280 to the original ERA5-Land resolution (0.1°; 11km). We then compare the results of the kriging 281 against two other commonly used interpolation techniques: nearest neighbour and bilinear 282 interpolation. The results are summarised in Figure 3. Figure 3A shows a box plot of the mean 283 absolute difference between the downscaled SAT and the original (i.e. upscaled) data. Nearest 284 neighbour interpolation is the simplest approach and so gives us a useful point of reference. 285 We see nearest-neighbour interpolation produces the largest differences, with bilinear 286 interpolation giving us slightly better results. The best interpolation is the kriging using just 287 elevation as a co-variate. When we start to add other co-variates, the estimates get further from This was combined with the dynamical uncertainty to create a total uncertainty, .We 293 compared this total uncertainty to the error in our estimated SAT by fitting a normal distribution 294 to the differences between the kriged SAT and the original ERA5-Land data for all months in 295 the period 1981-2010 to determine the spread in the downscaling error ( ). Figure 3B shows 296 the normalized difference between these two uncertainties. In all but a few coastal locations 297 the uncertainty given by KrigR is larger than the actual spread in the errors; meaning that in 298 almost all cases the difference between our downscaled data and the original ERA5-Land data 299 lies within the uncertainty from KrigR. 300 The picture is a little different for Qsoil. In Figure 3C we can see that a simple nearest-301 neighbour interpolation actually gives us the best results. Furthermore, kriging with any of the 302 covariates we considered actually produces worse results than a simple bilinear interpolation. 303 However, the big advantage with kriging is that we also get an uncertainty associated with our 304 estimate, and from Figure 3D we can see that the difference between the downscaling with 305 kriging and the original ERA5-Land data is always within the uncertainty given by KrigR. So 306 even in the case we do not have useful co-variates for our chosen downscaling, we can be 307 confident that the uncertainty from the kriging captures the real uncertainty in our downscaling 308 i.e., that the actual value at our target resolution, i.e. from the original ERA5-Land data, lies 309 within the statistical uncertainty given by KrigR.

TRADE-OFF WITH LOCAL KRIGING 320
One of the choices in KrigR is whether and how to localise the kriging. This defines the number 321 of neighbouring gridcells (nmax) used to derive the relationships between the field to be 322 downscaled and the co-variates. The larger the choice of nmax, the larger the number of cells 323 that will be used in the kriging process, and the closer the relationship between the target 324 variable and the co-variates becomes to the domain-average relationship. If nmax is set too 325 small then the relationships with the co-variates may be spurious, in which case so too will be 326 the downscaled product. This is represented in the fact that the smaller the nmax, the larger the 327 uncertainty in the downscaled product. This effectively enables users to fine-tune their choice 328 of nmax for individual needs and requirements within the cost-benefit trade-off between 329 computational resources and data uncertainty.