Long‐term rainfall regression surfaces for the Kruger National Park, South Africa: a spatio‐temporal review of patterns from 1981 to 2015

As an important bottom‐up driver of ecosystem processes, rainfall is intrinsically linked to the dynamics of vegetation and species distributions through its effects on soil moisture content and surface water availability. Rainfall effects are thus spatially and temporally specific to different environmental role‐players. Knowledge of its spatio‐temporal pattern is therefore essential to understanding natural ecosystem flux and potential climate change effects. Climate change poses a serious threat to protected areas in particular, as they are often isolated in fragmented landscapes and confined within hard park boundaries. In consequence, a species' natural movement response to resulting climate‐induced niche shifts is often obstructed. Long‐term, accurate and consistent climate monitoring data are therefore important resources for managers in large protected areas like the Kruger National Park (Kruger). In this article we model local rainfall measurements as a function of global rainfall surfaces, elevation and distance to the ocean using a generalized additive mixed effects model to produce fine‐scale (1 km2) monthly rainfall surfaces from July 1981 to June 2015. Results show a clear seasonal cycle nested within an oscillating multi‐decadal trend. Most noticeably, seasonality is shifting both temporally and spatially as rainfall moves outside of the typical dry/wet periods and areas. In addition, high‐rainfall seasons are generally receiving more rainfall while low‐rainfall seasons are receiving less. Northwestern regions of the park are experiencing more extreme annual rainfall differences, while far northern and southern regions show greater seasonality changes. The well‐described north–south and east–west rainfall gradient is still visible but the spatial complexity of this pattern is more pronounced than expected. Taken together, we show that Kruger's spatio‐temporal rainfall patterns are changing significantly in the short to medium term. The resulting raster data set is made freely available to promote holistic ecosystem studies and support longer‐term climate change research (http://dataknp.sanparks.org/sanparks/).


Introduction
Rainfall plays a central role in a myriad of natural processes, including river health, the transportation of nutrients (Strauch, 2013), soil moisture (Berry and Kulmatiski, 2017), vegetation dynamics (Ekblom et al., 2012), fire regimes (Archibald et al., 2009), animal movement and distribution patterns (Seydack et al., 2012) and landscape heterogeneity (MacFadyen et al., 2016). Within protected areas these processes function together to safeguard ecosystem integrity, helping to ensure the survival of places like the Kruger National Park (Kruger) as strongholds of biological diversity, socio-economic wealth and protectorates of natural resources and other ecosystem services (Secretariat of the Convention on Biological Diversity (SCBD), 2008). The spatio-temporal patterns of rainfall are thus an important variable to include in any ecological study, especially in the face of current climate change predictions (Hitz and Smith, 2004;Chapin III et al., 2011). The Intergovernmental Panel on Climate Change's (IPCC's) recent projections warn of severe declines in renewable surface and ground water resources for the subtropics (IPCC, 2014). South Africa (SA) in particular is projected to become progressively hotter and drier (Department of Environmental Affairs (DEA), 2013), raising serious concerns for local protected area management agencies to respond (Van Wilgen et al., 2016). These concerns stem from protected areas being especially vulnerable to the effects of climate change because they are often isolated in fragmented or bounded landscapes where species movement response is restricted (Hannah et al., 2007). Subtle changes to species (floral and/or faunal) RAINFALL PATTERNS OF THE KRUGER NATIONAL PARK 2507 composition and distribution dynamics in these natural systems can therefore act as early warning indicators of climate change effects (Garcia et al., 2014).
Ecosystem level research, practiced at the scale of a protected area, can then be an important addition to global climate monitoring efforts if channelled through long-term ecological research networks (Knapp et al., 2012). This type of research is often focussed around elucidating cause-and-effect relationships among various components of the ecosystem to anticipate potential state changes that may compromise system integrity (Pienaar, 1977). To do this effectively, reliable data on bottom-up substrates (e.g. topography, climate, geology, soils and vegetation) and top-down agents (e.g. fire and herbivory) are needed to discern undesirable change from natural ecosystem dynamics (Pickett et al., 2003;Venter et al., 2003). The need for robust, well-documented and continuous long-term climate data is therefore crucial (Kemp et al., 2012). Long-term data sets make it possible, for example, to quantify the effects of past extreme climatic events on ecosystem dynamics (e.g. vegetation patterns) to provide insights into the potential effects of varying climatic conditions in the future.
A large number of global climatological data sets have become available in recent years and are increasingly easy to access and use (see examples in Table S1, Supporting information). They are however modelled predictions of real-world climate, inherent with statistical bias and uncertainty (Muñoz et al., 2011;Kearney et al., 2014;Stoklosa et al., 2015). Moreover, the coarse resolutions at which these data are generated (e.g. 5 km CHIRPS; Table S1) are often not suitable for finer-scale ecological research (Kearney et al., 2014). To overcome these limitations a number of bias correction (Haerter et al., 2015) and spatial interpolation methods have been developed (Wagner et al., 2012). Most common among these, especially for data sparse regions like Africa, is the use of ground-based measurements and environmental covariates to augment globally modelled climate surfaces (Plouffe et al., 2015;Verdin et al., 2016). The importance of local long-term climate networks therefore needs to be more firmly recognized (Van Wilgen et al., 2016). Because the effects of rainfall on ecosystem dynamics play out at different spatio-temporal scales, local finer-scale data will remain irreplaceable sources of information for both regional-and global-scale ecosystems research Odiyo et al., 2015;MacFadyen et al., 2016). We illustrate this by providing a comprehensive analysis of regional (∼20 000 km 2 ), short to medium-term patterns and trends of local rainfall in Kruger from 1981Kruger from to 2015 Although countless studies conducted in Kruger cite the importance of climate, specifically rainfall, many are restricted to using park-wide averages or nearest station estimates (Owen-Smith and Ogutu, 2003;Smit et al., 2013;Trollope et al., 2014). We expect the dynamics of Kruger's rainfall are more complex, varying at different spatio-temporal scales that are undetectable when using simple regional averages. As a consequence, the spatio-temporal range of potential rainfall effects could also be more dynamic and further reaching than previously studied. We try to (1) answer some of these questions by assessing the temporal trends and spatial patterns of rainfall in Kruger over the last 34 years, (2) identify significant changes to the spatial rainfall patterns and ultimately (3) make the resulting spatio-temporal regression surface products freely available for future research.

Study area
Some of the world's largest protected areas are found in African savannas (IUCN and UNEP-WCMC, 2016). Large protected areas in these savannas, like Kruger, are important for the preservation of large-scale intact environments and the conservation of diverse wildlife assemblages and extensive wilderness qualities (Cantú-Salazar and Gaston, 2010). SA has a highly variable climate and pronounced rainfall seasonality, driven predominantly by sea surface temperatures, the effects of the El Niño-Southern Oscillation (ENSO) -both its warm (El Niño) and cool (La Niña) phases and the Intertropical Convergence Zone (ITCZ) (Philippon et al., 2014;Dedekind et al., 2016;Favre et al., 2016). The displacement of ITCZ to the south of the equator during the austral summer produces the South Indian Convergence Zone which restricts most of SA to a summer rainfall region (Dedekind et al., 2016). Kruger falls within a semi-arid zone which spans much of the central and northeastern parts of SA (Rutherford and Westfall, 1986;Trabucco and Zomer, 2009). On average, 84% of Kruger's total annual rainfall is concentrated between the summer months of November and April (Zambatis, 2003). Regionally, the area is split into two distinct climatic zones: Northern Arid Bushveld with 300-500 mm annual rainfall and Lowveld Bushveld with 500-700 mm annual rainfall, roughly north and south of the Olifants River, respectively ( Figure 1; South African Weather Services (SAWS), 1986). Locally, the terrain is moderately undulating (100-500 m a.s.l) with three higher-lying areas, namely the Shitshova range near Punda Maria in the northwest (650 m a.s.l), Lebombo mountain range bordering Mozambique in the east (480 m a.s.l) and Khandizwe, near Malelane in the southwest (840 m a.s.l; Figure 1). Here updrafts associated with higher elevation increase rainfall and with it soil leaching, thereby producing sandier more nutrient poor soils. For the rest of Kruger, rainfall is generally accepted to increase from north to south and from east to west as altitudes increase closer to the Drakensburg escarpment (see Venter et al. (2003) for a full account of local climate).

Local climate station records
Kruger has a long history of weather and rainfall monitoring, started by early naturalists in 1903 (Zambatis, 2003). Park rangers were later responsible for recording in situ rainfall measurements from 1920, which were augmented with six formal weather stations erected in partnership with the South African Weather Services (SAWS) in 1961 (Joubert, 2007). These stations included Skukuza, Pretoriuskop, Satara, Letaba, Shingwedzi and Punda Maria ( Figure 1). More stations were gradually added, culminating in 43 recording stations located across the park and providing a good representative sample of Kruger's different climatic zones, habitats and geographic range ( Figure 1). Apart from a number of sanctioned station closures over the years, missing data records threaten to compromise the quality and reliability of Kruger's climate station data. From the original 43 stations in Kruger we filtered for those stations with <10% missing daily rainfall records, resulting in 23 stations with a reliable data collection history from 1981 to 2015 ( Figure 1). We obtained additional rainfall records for 39 stations outside of Kruger, within 100 km of the boundary, with a similar rate of <10% missing daily records (South African Weather Services (SAWS), 2012). Another 12 stations from Kruger and 47 from SAWS, containing >10% but <20% missing values, were combined for cross-validation. In total, we collated daily rainfall records from July 1981 to June 2015 for 62 rain gauges for model training and 59 rain gauges for model validation (total n = 121). Rainfall records were reclassified into rainfall years (1 July-30 June) rather than calendar years to capture a full dry and wet season per rain year, i.e. in Kruger rainfall is generally restricted to December, January and February months, with little falling in June or July ( Figure S1).

Covariates: CHIRPS, DEM and DIOC
Daily and monthly CHIRPS GeoTIFF data were downloaded for the period July 1981 to June 2015 from the Climate Hazards Group (GHCN) data portal (Funk et al., 2015). CHIRPS is a gridded global rainfall surface product, which combines local station data with satellite rainfall estimates to produce a moderate resolution (∼5 km), long-term precipitation data set (Funk et al., 2015). Being a global product, Kruger's rain gauges are poorly represented in the GHCN's network of climate stations, i.e. only three stations inside Kruger's boundary and another two SAWS stations within 2 km of the park boundary are represented in GHCN (Menne et al., 2012). Resulting daily (n = 12 418) and monthly (n = 408) CHIRPS Geo-TIFFs were then stacked and the CHIRPS rainfall values extracted for each climate station locality (Hijmans, 2015).
Ground elevation values were also extracted from a 5 m digital elevation model (DEM) (Van Niekerk, 2012) for the same 62 sites and the straight line distance to the Indian Ocean coastline (DIOC) calculated. DEM and DIOC were included here because the combined effects of elevation and distance to the warm Agulhas current are expected to increase rainfall. More specifically, as the warm inland air meets the cooler sea air and the orographic effect of topography forces this air upwards, moisture condenses forming rain (Aalto et al., 2013).

General rainfall trends
Mean annual rainfall (MAR) summaries from 1910 to 2015 were extracted from Kruger's rainfall records, from which the long-term annual mean (541 mm) was calculated, along with a 3-year moving average. These summaries were derived from a variable number of stations across Kruger for different years; for example, MAR in 1910 was calculated using one station while 2015 comprised 22 stations. A continuous Morlet wavelet transform of mean monthly rainfall from July 1981 to June 2015 was then preformed to identify temporal frequency and variance patterns in Kruger's rainfall (Chan, 2000). Wavelet analysis is an established tool for climate research because of its ability to decompose non-stationary signals, like rainfall, into time-frequency representations using the wavelet power spectrum (de Jongh et al., 2006;Unal et al., 2012). The spectrum depicts the relative strengths (variance on z-axis) of rainfall frequencies (months on y-axis) smoothed across time steps (years on x-axis) to deliver a contour image illustrating rainfall periodicities (Chan, 2000).

Spatio-temporal regression surfaces
Missing daily rainfall records were first filled using a nonparametric imputation procedure with random forests (Stekhoven, 2013). In this way missing values from our 62 stations were filled using a random forest model with the observed daily rainfall measurements and related CHIRPS records (see Stekhoven and Buehlmann (2012) for full methodological detail). This process yields an 'out-of-bag' error estimate, for which our data had a 0.001% error rate expressed as a normalized root-mean-square error (NRMSE) of rainfall in millimetres. The filled data series was subsequently summed by year and month (year-month) to produce monthly summaries (n = 12) for each station (n = 62) and year (n = 34). A full 1 km 2 spatio-temporal grid (raster) was then generated over Kruger for 12 months by 34 years (n = 408) and the associated monthly CHIRPS, DEM and DIOC values extracted. Our response variable (rainfall) was therefore known at some dispersed localities (n = 62) across the study area for the study period (n = 408) while our explanatory variables (CHIRPS, DEM and DIOC) were known over the full spatio-temporal domain (Kilibarda et al., 2014). A correlation matrix confirmed no serious problems of collinearity between covariates that could cause overinflated standard errors ( Figure S2). Importantly, however, CHIRPS, DEM and DIOC are all spatial covariates and therefore may be collinear because of their geographic position, thus making them inherently spatially autocorrelated (see latitude and longitude entries in Figure S2). We accounted for both spatial and temporal autocorrelation by including a nested rational quadratic correlation structure (corRatio) which captured the dependence between observations close in space (u) and nested in time (t) (Hefley et al., 2017). Rainfall and thus its relationship with our covariates are also expected to be non-stationary. We accounted for this non-stationarity using a general additive mixed effects model (GAMM) which allows the relationship of our measured rainfall and CHIRPS to change smoothly across continuous levels of DEM and DIOC (spatially) and time (seasonally and annually). In this way our measured rainfall response is explained by its nonlinear relationship with CHIRPS, DEM and DIOC over space, seasonal cycle (within year) and longer-term trend (between years), including both fixed and random effects (Zuur et al., 2009).
Specifically, we used different cubic regression spline (CRS) smoothing bases to estimate rain (rainfall in mm) at location u (u = 1 to 62) and time t (t = 1 to 408). First, a standard CRS (f 1 ) was used for chirps (CHIRPS), dem (DEM) and dioc (DIOC). A cyclic CRS (f 2 ) was then used for mnth (month) with 12 pre-assigned knots, and again a standard CRS (f 3 ) for year (rainfall year) with 34 pre-assigned knots. The random effects caused by interactions between variables were estimated using tensor product interaction terms, namely (f 4 ) for year and mnth with CRS and cyclic CRS smoothers, respectively; (f 5 ) for chirps, dem and dioc with Gaussian process (kriging model) smoothers which account for the nonlinear effect DEM and DIOC have on the contribution of CHIRPS to the fitted rainfall response; and finally (f 6 ) for chirps, year and mnth with a thin plate regression spline, cyclic CRS and CRS smoothers to account for seasonal cycles and longer-term trends. represents i.i.d random variables expressed as N(0, 2 Λ) with Λ describing the residual autocorrelation nested within time t and location u using corRatio (Wood, 2006). A variety of model structures were tested until ultimately the model with the lowest Akaike information criterion (AIC) was selected as the best predictor of Kruger's spatio-temporal rainfall patterns (Zuur et al., 2009). Described above (Equation (1)), predictions generated from this model were based on the full 1 km 2 spatio-temporal grid of covariates for Kruger. Cross-validation of results was performed using a Pearson's product-moment correlation of predicted monthly rainfall surface values with ground measurements from 59 stations held back for model validation. The root-mean-square error (RMSE) is also reported, providing a good overall measure of model performance (Aalto et al., 2013). Final outputs were saved as 408 interpolated 1 km 2 surface grids describing Kruger's monthly rainfall patterns from 1981 to 2015. All analyses were conducted using statistical software R, version 3.1.3 R Core Team, 2016), and associated R packages referenced in the text.

Spatio-temporal variations in rainfall
We described the seasonal (within year) and longer-term (between years) trends in the resulting rainfall surface layers using another GAMM with spline interactions. Specifically, the within year seasonal effect (t 1 = 1 to 12 months) of our predicted monthly rainfall was allowed to vary smoothly with a between year trend effect (t 2 = 1 to 34 years) over space (u) (Wood, 2006;Zuur et al., 2009;Van Rij et al., 2015) as follows: √r whererain is the GAMM (Equation (1)) estimated measurement of rainfall at location u (latitude y and longitude x) and time t. The seasonal cycle mnth and overall trend year were modelled as f 1 , a cyclic CRS and f 2 , a standard CRS with 12 and 34 knots. f 3 , the tensor interaction thereof was included to reflect how rainfall responds to different combinations of month and year in comparison to the overall trend and seasonal cycle in f 1 and f 2 (Wood, 2006).
A thin plate CRS f 4 , of the x, y coordinates at location u, was also included to examine the distribution of monthly rainfall across Kruger's landscape. The amount by which this pattern differs over time was analysed using f 5 , a tensor interaction of x, y coordinates at location u with mnth and year at time t.

General rainfall trends
As expected, Kruger's rainfall patterns are cyclical in nature, varying seasonally from winter to summer and oscillating between wet and dry cycles ( Figure 2). Upon inspection of annual rainfall deviations from long-term means for the past 105 years, wet and dry cycles appear to occur about every 5 years (Figure 2(a)). From Figure 2(a) it appears these high-and low-rainfall periods generally match La Niña and El Niño years (see Figure S3 for Kruger's MAR in relation to ENSO events). Patterns become clearer with a continuous Morlet wavelet transform of mean monthly rainfall from July 1981 to June 2015 (Figure 2(b)). Results show the most significant pattern persists at an inter-annual scale of 10-16 months (Figure 2(b)). Short-term noise dominates monthly scale patterns with high-rainfall years blurring wet and dry season extents. For example, the extreme rainfall year 1999/2000, during which Kruger experienced severe flooding, shows significantly concentrated power (darker shades) regions already from the third month through to the twentieth (Figure 2(b)). With the occurrence of less extreme rainfall events over the last 10 years, inter-annual rainfall patterns are clearer. However, the power band between 10 and 16 months does appear to be widening slightly, suggesting seasonality may be shifting. A significant 4-5 year cycle (48-60 months) persists from   1986/1987/1988/1991/1992/1995/1996/1999/2000/2005 year cycle is however outside of the wavelet plot's 'cone of influence' and therefore requires more data to interpret correctly (Figure 2(b); Bunn et al., 2016). The spline smoothed trend line suggests no clear long-term trend is visible from 1981 to 2015 (Figure 2(b)).

Spatio-temporal regression surfaces
Model selection was based on AIC ranks, where the model with the lowest AIC was considered the best model of Kruger's rainfall (see Equation (1) for selected GAMM equation and Figure S4 for model diagnostics). The nonlinear trend of time and space with CHIRPS, DEM and DIOC and a temporally nested spatial autocorrelation structure explained 73% of the deviance (R 2 ) in Kruger's monthly rainfall (Table S2). The contribution of CHIRPS, DEM and DIOC was significantly nonlinear (p < 0.001) with all estimated degrees of freedom (edf) > 8.0, thus precluding the use of standard linear regression techniques (Zuur et al., 2009; Table S2). The smoothed seasonal and longer-term trend effects were also significant (p < 0.001; Table S2). CHIRPS was strongly positively related to measured rainfall (Figure 3(a)), while the effects of DEM were more significant at higher altitudes falling outside of Kruger's altitudinal range of 200-840 m (Figure 3(d)). DIOC elicited an unclear effect on rainfall (Figure 3(e)) but its interaction effect on CHIRPS' contribution to rainfall response was significant (p = 0.005; Figure 3(h)), as was that of DEM (p < 0.001; Figure 3(g)). In other words, the relationship of CHIRPS with station measured rainfall varied significantly as a function of elevation and distance to the coastline. A strong seasonal cycle was clearly visible (Figure 3(b)) with a more variable longer-term trend showing decadal fluctuations in rainfall response (Figure 3(c)). The nonlinear effects of CHIRPS also varied significantly with season (p = 0.025; Figure 3(f)) and year (p < 0.001; Figure 3(i)). Predicted rainfall surfaces, using GAMM (Equation (1)) and the full 1 km 2 spatio-temporal grid of associated covariates, showed visible spatio-temporal changes to seasonal and annual rainfall (Figure 4; see Video S1 for full animated results). Cross-validation of these regression surfaces against 59 model training stations (n = 7645) confirmed predictions were accurate overall (R 2 = 0.78, t = 107.5, df = 7643, p < 0.001) with an RMSE error of 19.99 mm year −1 . The R 2 and RMSE did however vary over time, with higher-rainfall years generally exhibiting smaller R 2 variations and greater RMSE variations ( Figure S5).

Spatio-temporal variations in rainfall
There exists a clear seasonal cycle (p < 0.001, edf = 10; Figure 5(a)) with strong evidence of longer-term patterns that vary significantly over 10 and 20 year cycles, making it difficult to interpret a general overall trend (p < 0.001, edf = 9; Figure 5(b)). The interaction of rainfall year and month was significant (p < 0.001; edf = 90) and showed patterns of extreme rainfall conditions (Figure 5(c)). In fact, taking a ruler to Figure 5(c), beginning between August (8) and September (9), one could almost draw a perfect diagonal of high-rainfall period moving later in the season until about April (4) 2007. Periods with earlier rainfall appear in sync with documented La Niña events while the El Niño link to lower-rainfall periods is less clear (Figures 5(c) and S3).
Rainfall response also differed significantly (p < 0.001, edf = 29) across space, with the effect of location (x, y) being significantly nonlinear (Table S3) as rainfall varied from north to south and west to east ( Figure 5(d)). Rainfall was highest in the far southwest and lowest in the mid-northwestern parts of Kruger ( Figure 5(d)). These patterns however changed annually and seasonally ( Figure 5(e)). Distinct periods of high rainfall were visibly more pronounced in the south compared to the north (Figure 5(e) and iii). Similarly, higher-rainfall periods were more pronounced in the west compared to those experienced in the east (Figure 5(e)-ii and iv; see Table  S3 for a full description of model results). These results suggest more localized or regional changes may dominate Kruger's rainfall dynamics than previously studied.

Discussion
A recent assessment by the IPCC rated Africa as having a 'high and significant' risk of climate change, and thus should expect severe water resource constraints in the near future (IPCC, 2014). This desiccation is further supported by Chen and Chen (2013) and Rubel et al. (2017) in their reassessments of Köppen (1930) global climate classification for the period 1901-2010, which shows much of Kruger becoming drier. Adding to these results, we found a more noticeable spatio-temporal change in the short to medium-term seasonal cycles of low-high rainfall, rather than a clear increasing or decreasing trend. Our findings are supported by a recent larger-scale assessment of climate change across all South African National Parks (SANParks) by Van Wilgen et al. (2016), who found increased rainfall variability in parks to the east. Similarly, they found no significant annual trend in Kruger's rainfall but suggest an increase in seasonality with longer dry periods (Van Wilgen et al., 2016). At a finer scale we found this shift to be spatio-temporally dynamic. Figure 5(c) illustrates the seasonal shifts in rainfall patterns as they are affected by strong La Niña events. Philippon et al. (2014) found similar patterns for southern Africa as vegetation greenness (NDVI) was dampened during El Niño and increased during La Niña events. Figure 5(d) shows the distribution of rainfall across the park from west (e.g. Pretoriuskop) to east (Mozambique border) and from south (e.g. Malelane) to north (e.g. Pafuri). As described by Gertenbach (1980) and Zambatis (2003), the more southern latitudes generally receive higher rainfall, as do the more western longitudes. However, Figure 5(d) suggests this spatial pattern is clearly more complex than can be adequately accounted for by regional averaging of the past. Moreover, these spatial patterns are also continually changing with both high-and low-rainfall periods lengthening or shortening, and/or shifting outside of recognized wet/dry season windows over the short to medium term ( Figure 5(e)). For example, the northern latitudes have gone through extreme wet, dry, wet cycles in comparison to the south from 1981 to 2015 ( Figure 5(e)). This pattern is mirrored in the western longitudes ( Figure 5(e)ii), suggesting the far northwest of Kruger is undergoing significant spatio-temporal rainfall changes over the short to medium term.
A strong seasonal shift is also apparent in the west, where October (early season) rainfall is generally increasing while June (mid-dry season) is decreasing ( Figure 5(e)iv). Similarly, in the far north and south of Kruger, summer rainfall is generally increasing while winter rainfall is decreasing (Figure 5(e)-iii). In contrast, Kruger and Nxumalo (2017) analysed the historical rainfall trends in SA from 1921 to 2015 and found no significant trends in seasonal rainfall totals for our area. They did however find a significant downwards trend in total annual rainfall for the northern and southern regions of Kruger, although this signal was not significant for individual stations in the area (Kruger and Nxumalo, 2017). Our results suggest that these linear averages conceal seasonal and spatial rainfall trends at the local level in Kruger. If wetter seasons and areas are becoming wetter while drier seasons and areas are becoming drier, as seen in Figures 5(c)i-iv, park-wide averages would indeed obscure these changes. Taken together these shifts in extremes could lead to an increased risk of droughts or floods (Pohl et al., 2017), like  (b) long-term trend in rainfall with oscillating high-and low-rainfall periods; (c) contour plot of the partial nonlinear interactions of rainfall year and month showing the amount by which the fitted monthly rainfall is adjusted from the overall trend and seasonal cycle for any combination of month and year. Strong El Niño years are marked with _ below the year label and strong La Niña years withabove year on the x-axis; (d) perspective plot of the interaction effect of longitude (x) and latitude (y) on rainfall response; (e) partial contour plots of latitude (south-north) and longitude (west-east) tensor product interaction terms showing the amount by which the fitted monthly rainfall is adjusted from the overall trend and seasonal cycle from south to north (i, iii) and west to east (ii, iv). [Colour figure can be viewed at wileyonlinelibrary.com].
the floods experienced in Kruger in 2000(Knight and Evans, 2017. Although these spatio-temporal changes may obscure broad-scale trends (e.g. Kruger and Nxumalo, 2017;Van Wilgen et al., 2016), they will continue to have significant effects on ecosystem processes like species population dynamics and fire. As a result, they should form an important component for management planning. For example, in understanding these spatio-temporal dynamics, park management can better explore the effects of different climate change forecast scenarios, by examining the environmental response to past extreme climatic events. The severe drought of 1991/1992 for example had significant effects on woody and herbaceous plant species composition and consequently predator-prey dynamics, as well as a negative effect on river geomorphic diversity (Heritage and van Niekerk, 1995;Viljoen, 1995;Zambatis and Biggs, 1995). Similarly, the large Sabie River flood in February 2000 changed river geomorphology and associated riparian vegetation and species distribution patterns (Parsons et al., 2006). Future studies might consider including the new high-resolution (1 km 2 ) global cloud cover product described by Wilson and Jetz (2016). This new data set combines moderate resolution imaging spectroradiometer (MODIS) imagery from 2000 to 2014 to derive long-term cloud-cover dynamics aimed at improving spatial autocorrelation errors, typical of interpolated climate data, towards a better spatio-temporal understanding of global biodiversity and ecosystem properties (Wilson and Jetz, 2016).
Rainfall is indisputably a significant driver of numerous bottom-up and top-down drivers of biodiversity in protected areas (Garcia et al., 2014). It is strongly associated with vegetation dynamics and thus ungulate population dynamics (Owen-Smith and Ogutu, 2003;Dunham et al., 2004), predator densities and demography (Celesia et al., 2010) and disease outbreaks (e.g. malaria, Colón-González et al., 2016). Reductions in rainfall can have severe implications for the provision and management of artificial water, which in itself may prove to have far-reaching impacts (Smit, 2013). In light of the landscape-scale changes we found, research and monitoring in Kruger will need to become more spatio-temporally dynamic with climate considered more as a multi-scale driver than a general stratification template.

Conclusion
In Kruger, the importance of rainfall monitoring was stressed as far back as 1901, which is reflected in the current research objectives describing rainfall as a major factor effecting biodiversity (SANParks, 2008). Despite its obvious importance, funds for inventorying and long-term monitoring are increasingly restricted as focus shifts to other conservation priorities (e.g. anti-poaching). Amplified by the 'publish-or-perish' mentality and the need for rapid data collection and analysis, the time-cost-benefit of maintaining long-term in situ monitoring programs continues to wax and wane (Sarewitz, 2016). However as we have shown, local rain gauge observations in combination with satellite-derived rainfall estimates can both supplement global products for local finer-scale studies as well as potentially improve the accuracies of global assessments. To support these initiatives, we provide our data as an open source gridded rainfall surface data set for Kruger from 1981 to 2015 to encourage the development of a better understanding of rainfall-driven system response patterns (http://dataknp.sanparks.org/sanparks/).
In light of projected climate change impacts, protected areas in particular may become increasingly susceptible as hard park boundaries act as physical barriers to animals seeking to move in response to bioclimatic range shifts. Without adaptive conservation strategies and climate-change conscious spatial planning, these systems may ultimately transform into undesirable states as a consequence of altered species dynamics. Because the effects of altered climatic conditions will be landscape and species specific, studies of local climate change patterns need to be spatio-temporally explicit. Rainfall, for example, will act as a driver of change at multiple scales and should be assessed as such when used to make management decisions. Parks for access to rainfall station data. This project was supported by the National Research Foundation (NRF) of South Africa (grant nos. 89967, 109244 and 109683). Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF.

Supporting information
The following supporting information is available as part of the online article: Table S1. Examples of some global climatological data sets that have become available in recent years (Muñoz et al., 2011;Kearney et al., 2014See NCAR (2014 and Sun et al. (2016) for a full summary. Table S2. GAMM results describing the response of ground measured rainfall as explained by its nonlinear relationship with CHIRPS, elevation (DEM) and distance to the Indian Ocean coastline (DIOC) smoothed over space, seasonal cycle (within year) and longer-term trend (between years). Table S3. GAMM results describing the seasonal (within year) and longer-term (between years) trends in predicted rainfall surface layers, allowing the within year seasonal effect (mnthRain) to vary smoothly with a between year trend effect (yearRain) over space (x, y). Figure S1. Level plot of month in which the maximum daily rainfall was received for a particular rainfall year and selected stations in Kruger from 1981 to 2015. The maximum value is displayed in each year-month-station cell. Figure S2. Visualization of correlation matrix confirming no problems of collinearity between variables but highlighting inherent spatial autocorrelation between spatial covariates. Figure S3. The mean annual rainfall (MAR) of the Kruger National Park and the Oceanic Niño Index (ONI) identifying El Niño (warm) and La Niña (cool) events in the tropical Pacific from 1950 to 2015. Years labelled red represent strong to very strong El Niño events, which typically decrease rainfall across southern Africa. Years labelled blue indicate strong La Niña events, where rainfall is generally increased (Huang et al., 2017 Data from http://origin.cpc.ncep.noaa.gov/products/analysis_ monitoring/ensostuff/ONI_v5.php (accessed 19 October 2017). Figure S4. Model diagnostics for GAMM (Table S2; see Equation (1) for full GAMM equation) showing no residual temporal autocorrelation. Figure S5. Cross-validation of gridded results with monthly rainfall data from 59 stations (n = 7645) confirmed our spatio-temporal regression surfaces produced accurate predictions of monthly rainfall (R 2 = 0.78, t = 107.5, df = 7643, p < 0.001) with an RMSE error of 19.99 mm year −1 . Video S1. Animation of predicted rainfall surfaces, using the results of GAMM and the full 1 km 2 spatio-temporal grid of associated covariates, illustrating the seasonal and annual spatio-temporal rainfall dynamics in Kruger from 1981Kruger from to 2015