Geoscientific Model Development TopoSUB : a tool for efficient large area numerical modelling in complex topography at sub-grid scales

Mountain regions are highly sensitive to global climate change. However, large scale assessments of mountain environments remain problematic due to the high resolution required of model grids to capture strong lateral variability. To alleviate this, tools are required to bridge the scale gap between gridded climate datasets (climate models and reanalyses) and mountain topography. We address this problem with a sub-grid method. It relies on sampling the most important aspects of land surface heterogeneity through a lumped scheme, allowing for the application of numerical land surface models (LSMs) over large areas in mountain regions or other heterogeneous environments. This is achieved by including the effect of mountain topography on these processes at the sub-grid scale using a multidimensional informed sampling procedure together with a 1-D lumped model that can be driven by gridded climate datasets. This paper provides a description of this sub-grid scheme, TopoSUB, and assesses its performance against a distributed model. We demonstrate the ability of TopoSUB to approximate results simulated by a distributed numerical LSM at around 10 4 less computations. These significant gains in computing resources allow for: (1) numerical modelling of processes at fine grid resolutions over large areas; (2) efficient statistical descriptions of sub-grid behaviour; (3) a “sub-grid aware” aggregation of simulated variables to coarse grids; and (4) freeing of resources for computationally intensive tasks, e.g., the treatment of uncertainty in the modelling process.


Introduction
Mountain regions extend over a large portion of the global land area and significantly influence climate as well as human livelihoods (Barnett et al., 2005;Gruber, 2012;Immerzeel et al., 2010).Complex topography in mountain regions causes high lateral variability of the surfaceatmosphere boundary by: (a) altering the local energy and mass fluxes between the ground and the atmosphere (caused by e.g., air temperature, shading, precipitation gradients); and (b) influencing subsurface properties (e.g., exposed bedrock in steep slopes, fine sediments and abundant water in valleys).Mountain environments are currently undergoing rapid and significant change worldwide due to global changes in the earth's climate e.g., warming mountain permafrost (Harris et al., 2003;Isaksen et al., 2001); retreat of mountain glaciers (Paul et al., 2007;Zemp et al., 2006;Barry, 2006); and reduction of snow cover in many regions (Laternser and Schneebeli, 2003;Mote et al., 2005).In order to understand the impact of these changes, under current and future climate conditions, tools are required to enable numerical modelling of physical processes that occur across a range of spatial scales.
Global climate models (GCMs) and regional climate models (RCMs) are able to generate continuous and physically consistent fields of climate variables for the observational period and for scenarios of future climatic conditions.However, the coarse grids such models operate on (∼ 10-500 km) limit the ability to resolve the interactions and feedbacks between the land surface and climate systems in complex topography, which is characterised by strong variability and nonlinear processes at the sub-grid scale (Giorgi and Avissar, Fig. 1.The scale problem TopoSUB addresses -(A) the complex geometry of mountain topography is aggregated to a mean value in coarse grids which does not necessarily account for sub-grid heterogeneity (B) application of numerical models requires a fine grid in order to account for the effect of this heterogeneity, which is computationally expensive.TopoSUB allows for application of numerical models over large areas through a lumped approach that samples the most important aspects of this heterogeneity.1997).A wealth of surface models exist which are capable of simulating processes in mountain regions on fine grids (∼ 1-100 m) and can be driven by coarse grid data with suitable regionalisation techniques (e.g.Bartelt and Lehning, 2002;Gruber et al., 2004;Klok and Oerlemans, 2002;Paul and Kotlarski, 2010), however, the strong fine-scale variability of mountain systems (cf.Gubler et al., 2011;Riseborough et al., 2008) precludes the application of high-resolution, distributed numerical models, over large areas.Therefore, the problem remains that despite the near-global availability of high-resolution digital elevation models, global climate datasets and suitable numerical simulators, land surface processes in complex topography remain poorly quantified in many aspects (Fig. 1).
This problem of scale has been previously approached through various forms of sub-grid parameterisation.This term can be defined as capturing the spatial variability of a modelled process at a suitable resolution, while reducing the demands for data and computation, by approximating its fine-scale distribution at a lower resolution (e.g.Hebeler and Purves, 2008;Giorgi and Avissar, 1997).Previous sub-grid approaches can be broadly classified as either discrete mosaic types or continuous probability density function (PDF) schemes (Wood et al., 1988;Avissar, 1991;Giorgi and Avissar, 1997).This distinction can also be conceptualised as the modelling of sub-grid inter-patch and intra-patch heterogeneity (Giorgi and Avissar, 1997).In mosaic approaches, a number of homogeneous subregions ("tiles") are defined at the sub-grid scale, each with its own energy, momentum and water budget.The surface fluxes are computed separately for each tile.Aggregation to the coarse grid is performed by averaging over the tiles which are weighted by their fractional cover (Avissar and Pielke, 1989;Koster and Suarez, 1992).Models differ on how these tiles are discretised.For example, Seth et al. (1994) and Dimri (2009) used a regularly spaced finescale sub-grid.Alternatively, a series of discrete classes based on surface vegetation type (Avissar and Pielke, 1989;Koster and Suarez, 1992) or topographical elevation (Leung and Ghan, 1995) have been used.Kotlarski (2007) developed a dynamic mountain glacier sub-grid parameterisation for inclusion in RCM's which explicitly accounts for run off generation and adjusts glacier area (dynamic tile) based on accumulation/ablation conditions.PDF-based approaches attempt to describe the variability of sub-grid characteristics through analytical or empirical distribution functions (Avissar, 1991;Famiglietti and Wood, 1994;Liang et al., 2006).This is based on the assumption that surface characteristics as well as climatic forcings vary according to distributions that can be approximated by the given PDF.This approach then explicitly calculates gridbox average surface fluxes for a surface variable distribution using numerical or analytical integration over the appropriate PDF.Walland and Simmonds (1996) used sub-grid statistics (variance, kurtosis) of distributions of topographical parameters to improve the simulation of the snowpack in GCMs.
Mountain regions exhibit more relevant dimensions which control land surface processes than more gently inclined areas (e.g., elevation, aspect, slope, etc.).This means that a simple mosaicing of the land surface is not appropriate and a more sophisticated technique is needed to account for this heterogeneity.While, statistical models exist (e.g.Boeckli et al., 2012) which provide good 2-D representation of single phenomena in complex topography, they usually do not resolve transient changes.Finally, methods exist (e.g., SAFRAN-Crocus scheme; Durand et al., 1993Durand et al., , 1999) ) which classify topography according to fixed classes based on terrain parameters and enable application of numerical models over large areas in a semi-distributed fashion.However, we suggest that a more flexible method that allows for continuous and unprescribed classification, which is optimised to the region and variables of interest, is a useful extension of such methods.
In this paper, we describe a sub-grid method which samples the most important aspects of land surface heterogeneity based on input predictors (PREDs) which describe important dimensions of variability in complex topography (e.g., elevation, aspect, slope, sky view factor).A lumped scheme then allows for the efficient application of numerical land surface models (LSMs) over large areas.Aggregation of simulated target variables (TVs, e.g., ground surface temperature, snow water equivalent) to the coarse grid and spatialisation to the fine grid is achieved through membership functions.We use cumulative distribution functions (CDF) of TVs to provide statistics of sub-grid behaviour over large areas.This Topographic SUBgrid tool (TopoSUB) allows for: (1) modelling of processes at fine grid resolutions, (2) efficient statistical descriptions of sub-grid behaviour, (3) a "sub-grid aware" aggregation of simulated TVs to coarse grids and (4) enables validation of results with fine-scale ground truth.The strength of the scheme is its ability to enable the computational representation of fine-scale processes over large areas by several orders of magnitude faster than a distributed model and, therefore, free resources for spatially or temporally expensive simulations as well as for exploring uncertainties in input data (e.g., climate projections) or LSM parameters and physics.While we acknowledge that a lumped approach compromises on 2-D representation (e.g., snow redistribution, surface runoff) it enables the application of sophisticated 1-D physics over large areas.
This paper provides a proof of concept of this tool by describing the method, providing guidance on parameter selection and performing validation experiments against baseline distributed model simulations.Whilst TopoSUB is designed as a tool for use in complex topography -the concept may also be of interest outside of mountain environments where alternative dimensions of variability are important.

K-means clustering
Samples are formed using the K-means clustering algorithm of Hartigan and Wong (1979), an unsupervised learning algorithm which is a suitable for clustering multidimensional data.K-means aims to partition all points into K clusters such that the total sum of squares (or squared deviations) from individual points (pixels) to the assigned cluster centroids in multivariate attribute space is minimised.This then represents an optimal clustering of points in attribute space for a prescribed number of samples.The algorithm is composed of the following steps: 1. Randomly place K points into the space represented by the objects that are being clustered.These points represent initial group centroids/seeds.2. Assign each object to the group that has the closest centroid.
3. When all objects have been assigned, recalculate the positions of the K centroids.
4. Repeat Steps 2 and 3 until the centroids no longer move.This produces a separation of the objects into K clusters from which the metric to be minimised can be calculated.
Although it can be proven that the procedure will always terminate, the K-means algorithm is sensitive to the initial configuration of cluster seeds and does not necessarily find the optimal configuration (Kanungo et al., 2002), corresponding to the global objective function minimum (as opposed local minima).The K-means algorithm is run multiple times to reduce this effect.
This multiple K-means algorithm has three important controlling parameters, K number of clusters, iter.maxmaximum allowed iterations of the algorithm and nstart number of random starts.Fixing K at 128 samples, a sensitivity analysis was performed on iter.max and nstart in order to define baseline parameter values to be used in TopoSUB.Stable performance, as measured by within sum of squares (WSS) was achieved for iter.max= 20 and nstart = 10 (Fig. 2).Iter.max shows a much greater variation in WSS over its range as it is an intrinsic part of the K-means algorithm.K-means is significantly sensitive to initialisation i.e., location of initial cluster seeds.For this reason it is highly recommended to run the K-means algorithm with several random starts and average the results.A first run of K-means algorithm is performed on a subset of input data (10 5 pixels) with 10 random starts and maximum iterations set to 20 (as previously defined).The cluster centres defined by this clustering of the subset are used to initialise K-means for the entire dataset (10 6 or more pixels).This allows for significant speed up of the algorithm (factor of 10, Table 1) while not compromising on the quality of results.

Fuzzy membership
In contrast to crisp membership (yes/no), fuzzy methods allow for varying degrees of membership to multiple sets (Zadeh, 1965).A membership function is given with a value in the interval 0-1.We apply fuzzy membership functions to allow for varying degrees of membership of pixels to multiple samples.Fuzzy methods of classification are usefully applied to multivariate classification problems when class overlap is required to represent continuous phenomena.The primary advantage of this method over crisp classification is that fuzzy methods allow for high resolution mapping of the TV, whilst accounting for topographic variability within each sample.Membership functions are calculated in two steps.First, the standardised squared distance (d 2 ) of the i-th pixel from the n-th sample centroid C of the j -th PRED is determined by: where Sd is class standard deviation (Burrough et al., 2001).
Then we can derive the membership, µ of the i-th pixel to the n-th sample using the formula of Sokal and Sneath (1967): where M is the so-called fuzzy exponent parameter.This parameter is a weighting exponent and it controls the degree of fuzziness of the membership grades.As M approaches 1, the clustering becomes crisper.As M becomes very large (i.e., M > 100), membership becomes almost constant so that clusters can no longer be distinguished.The determination of an optimal value for M in a fuzzy classification process remains an open question (e.g.Okeke and Karnieli, 2006;Burrough et al., 2000).Commonly cited values in the literature range from 1.3 to 3 depending upon application.Okeke and Karnieli (2006) propose a linear mixture model approach to optimise the value of M for any given dataset.We found an optimum value to exist between crisp (M = 1) and very fuzzy (M > 2) by iterating through values of M 1-2 in stages of 0.1.We identified 1.4 to give the most accurate results in majority TVs tested (compared with a distributed simulation), in our specific test case.However, the optimal value of M is not the focus of the present study and will not be further discussed.A standard fuzzy membership algorithm will compute memberships for all pixels to all samples, resulting in a membership matrix of n samples by p pixels.Even with the modest large area simulation we give in this paper, this results in matrix of magnitude 10 8 .In order to reduce storage demands we allow a reduced number of membership dimensions to be prescribed whilst preserving the functionality of the algorithm.This is achieved by either (a) prescribing absolute number of membership values allowed per pixel or (b) setting a target membership level (i.e., 95 %) and averaging the required number of memberships per pixel to achieve this threshold.

Statistical methods
Test distributions are compared with baseline distributions using a one-sample Kolmogorov-Smirnov test (KS-test) which is a non-parametric test for the equality of continuous, one-dimensional probability distributions that can be used to compare a sample with a reference probability distribution.If the sample comes from distribution F (x), then the KS statistic, D converges to zero.Comparison of simulated fine grid data against baseline grids is performed using the normalised root-mean-squared error (NRMSE) to allow comparison across simulated TVs of different scale ranges.The RMSE is normalised against the standard deviation of the observations, being a robust statistic that is less influenced by outliers than the range.Correlation statistics stated are Pearson product-moment correlation coefficients (r-value).

TopoSUB methods
TopoSUB is designed to provide an effective approximation of a spatially distributed grid with a lumped model.A key design principle is to be generic, allowing for choice of driving inputs (station data, gridded climate data), numerical model and output (TVs, resolution) to support a wide range of possible applications.It has two main modules (Fig. 3): a preprocessor to run only once and a post-processor to run many times together with an LSM.Because it is intended for use in  mountain areas, it needs the ability to accommodate more dimensions of variability than is usual in mosaicing techniques used to partition the sub-grid.Besides differing surface and subsurface properties, the effects of elevation, slope exposition, slope angle and horizon elevation are likely to be of importance.To allow the scalable use of this scheme, i.e., its application over large mountain ranges, it should employ repeatable and robust methods to form sub-grid samples.Because the influence of predictor variables (dimensions of subgrid variability that are accounted for) is not known a priori and may change laterally, a method for informed sampling is required in which the importance of predictor variables on the simulated target variable(s) is evaluated.Model parameters are given in Table 2 together with default values.

Pre-processor
The pre-processor configures the sub-grid by creating samples for later simulation in the LSM and by determining the membership functions of original 2-D pixels to those samples.For a given experimental domain this module need only be run once.The scheme can use several PRED variables from which the sub-grid scheme is constructed and usually these include DEM-derived land surface parameters.

Training routine
The input PREDs are initially clustered using the K-means algorithm to form a predefined number of samples.Scaling of PREDs is important because it affects the relative number of samples formed along a given dimension.With the scaling of PREDs, we, thus, influence how finely samples are resolved in which direction, a process that is important to optimise the number of samples.No scaling could result in, for example, differing results with elevation provided in units of metres or kilometres.At this initial stage, a simple scaling (centred and scaled, Eq. 3) is, thus, applied which normalises all PREDs to a standard scale under the assumption that all are of equal importance with respect to simulated TV.Sample centroids define the topographic and environmental input to the LSM that is run for this initial set of samples in a training simulation.

Sample formation by informed clustering
Based on the training routine results, one linear regression model is constructed for each of n TVs using the whole set of i PREDs as regressors (Eq.4) using generalised least squares (GLS).GLS is able to handle PREDs with non-normal distributions and/or which are partially correlated and, therefore, is more robust when implemented as an automated method.GLS minimises the squared Mahalanobis distance as opposed the residual sum of squares, as in regression methods using ordinary least squares.
The resulting regression coefficients, β i provide an informed scaling of the PRED i for the clustering algorithm by transforming them into equivalents of the TV n dimension and unit (Eq.5).As an example, elevation (m a.s.l.), slope angle ( • ) would be scaled into equivalents of ( • C) if ground temperature was the target variable. (5) Normalised to a sum of one, these coefficients can be averaged (β mean ) to accommodate more than one TV.Parameter W TV is an optional weighting function for individual TVs if extra information exists to justify higher weighting of a TV in the averaging process: TV = 1, ...n.
The set of PREDs, with informed scaling, are re-clustered to provide a predefined number of samples that are then more effectively distributed with respect to the desired TVs (Fig. 4).The sample centroids now provide the required input to the LSM.Additionally, the regression model can be optimised by using disaggregated r 2 (Genizi, 1993) either manually or automatically, by removing PREDs that contribute less than a stated threshold to model significance.
We have also implemented a routine which allows for the dropping of samples that are deemed insignificant.Insignificant samples are defined as those with fewer members than a percentage threshold (defined a priori) of a theoretical sample membership value obtained if all pixels where distributed equally among samples.The members of dropped samples are redistributed on a nearest neighbour basis to remaining samples, in euclidean space.

Pre-processor output
A vector of sample weights is calculated according to total membership of pixels to each sample.This provides the means by which to (a) aggregate TVs to the coarse grid and, (b) provide a rapid statistical description of sub-grid behaviour using a CDF.A matrix of membership functions of individual pixels to samples provides a means of spatialising results to the fine grid.We employ crisp and fuzzy membership.Crisp membership implies that pixels may belong to only one sample (that which they are assigned in clustering).Fuzzy membership allows for varying degrees of membership to multiple samples.This then accounts to some extent for within sample variance (in terms of PREDs) that inevitably exists and provides an alternative and smoother means of spatialisation at reasonable computational cost.
The final output from the pre-processor is the sub-grid configuration (Fig. 5) as defined by: (a) a small matrix of sample characteristics.These are the environmental characteristics of each sample used to drive the LSM as well as the aggregated sample weights used for aggregation to grid level and statistical description of sub-grid behaviour.(b) Membership information for the spatialisation of results to the fine grid.For crisp membership this has the dimensions of the original fine grid (B1), for fuzzy membership the ID and weight for the s most important samples are stored for each original pixel (Table 2), increasing the dimensions of this information to 2* s times the original size (B2), because the membership values as well as the sample ID's need to be stored per pixel.Based on the sub-grid configuration, the LSM is run in 1-D mode for each sample.

Post-processor
Based on the sub-grid configuration, the LSM output is postprocessed.The result of this can either be (a) summary statistics with respect to the coarse grid describing its sub-grid variability by a CDF or derived quantities; or (b) data spatialised to the original fine grid; or (c) data estimated for a list of discrete points to support validation studies using ground truth data.The coarse-grid summary statistics are computed according to the aggregated membership functions of individual pixels to each sample.TVs are spatialised to fine grid resolution according to the membership functions (crisp or fuzzy) of each pixel to each significant cluster.This accounts for the sub-grid heterogeneity that exists between cluster centroids.

Data and tools
The high-resolution input into the scheme is a 25 m digital elevation model (DEM, obtained from Swisstopo).The DEM used in this study covers a test area of 25.6 × 25.6 km (∼ 10 6 pixels) in the south-east of Switzerland (Fig. 6).The study area represents a good example of strongly variable mountain topography (elevation range 1556-4043 m a.s.l.) with which to test the performance of the scheme.Land

Land surface model
We employ the open-source LSM GEOtop (Endrizzi and Marsh, 2010;Dall'Amico et al., 2011;Rigon et al., 2006) which is a physically-based model that simulates the coupled energy and water balance with phase change in soil, a multilayer physically-based snow-pack model and surface energy fluxes in 1-D and distributed 2-D modes.It has been designed specifically for application in mountain regions.The model domain consists of a soil column of user-specified depth (typically of a few metres) from the ground surface, which is, in turn, defined by a Digital Elevation Model (DEM).The heat and subsurface water flow equations are then solved with finite differences schemes.
The multi-layer snow pack scheme accommodates compaction as well as water percolation and refreezing.The influence of topography on micro-climatology is parameterised, allowing for the solution of the surface energy balance for differing topographic situations based on one driving climate time series (Endrizzi and Marsh, 2010;Liston and Elder, 2006).A vegetation canopy was not considered in these experiments.The soil is uniform over the entire simulation domain, parameterised using the Van-Genuchten model (Van Genuchten, 1981), and has 5 layers and a total depth of 3.1 m.The model is run on an hourly timestep.We apply two years of spin up and then generate 1 yr of data.
The meteorological data, input as point time-series are spatially distributed by GEOtop to each simulation point using principles of the Micromet model (Liston and Elder, 2006).Specifically: (1) Air temperature follows a mean lapse rate (6.5 • C km −1 ).( 2) Since relative humidity is a nonlinear function of elevation, the dewpoint temperature is used for vertical extrapolation.(3) Model time step is used to calculate the solar radiation for that specific time.In addition, the influence of cloud cover, direct and diffuse solar radiation, and topographic slope and aspect on incoming solar radiation is accounted for.The distributed version has self and cast shadowing based on DEM, point has self and a uniform horizon elevation.(4) Precipitation is not adjusted.

Testing strategy
The primary aim of this evaluation is to test how well Topo-SUB is able to reproduce results of a distributed model.This is done by comparing results obtained from TopoSUB with results from a distributed LSM simulation on a regular 2-D grid.Both runs use the same LSM, GEOtop and same meteorology distribution scheme described in Sect.4.2.The key difference is that the simulation units are (a) samples resulting from clustering of predictors in TopoSUB and, (b) pixels at DEM resolution in distributed runs.The distributed runs have diverse spatial resolutions ranging from 25 m (10 6 cells) to 25 km (1 cell).TopoSUB is also run at differing resolution i.e., different levels of detail (number of samples) in the sampling of the sub-grid.In all experiments, both Topo-SUB and distributed runs are evaluated with respect to a 25 m distributed simulation (BASE), which we consider to be the baseline in this experiment, being the most finely discretised representation of the experiment domain.Experiments presented in results are: (Sect.5.1) grid aggregated results compared directly with corresponding statistics of BASE.For this, the mean and standard deviation as well as the 25th and 75th percentiles are calculated; (Sect.5.2) statistical description of sub-grid behaviour is evaluated by comparing the CDF of TopoSUB simulation to CDF of BASE simulations using KS-test; and, (Sect.5.3) fuzzy spatialised results are compared to BASE using the r-value and NRMSE to assess the predictive power of the scheme at the grid cell level.Finally, other aspects of TopoSUB performance are presented (Sects.5.4-5.6).To keep computation times for the 25 m resolution BASE simulation reasonable, water movement in the soil was not considered.

Model settings
For testing, input PREDs of elevation, aspect, slope and sky view factor are used in the clustering algorithm, all computed from the input DEM at 25 m resolution.In all experiments we simulate output response variables of air temperature (Tair, • C), ground temperature at 10 cm depth (GST, • C), snow water equivalent (SWE, mm) and incoming shortwave radiation (SWin, W m −2 ).These were chosen as suitable variables with which to test the performance of the TopoSUB scheme in terms of representing both near surface processes and energy fluxes.Air temperature represents a simple check of the scheme due to straightforward relationship with elevation in this study (using a standard lapse rate).SWin represents how well the scheme is able to represent topography with respect to parameters relevant to radiation modelling.Snow water equivalent and GST both represent important physical processes in mountain areas.TVs are analysed as mean annual values in all cases.Spatialised TopoSUB results are evaluated on a pixel by pixel basis whereas aggregated results are by definition evaluated as a mean value of the simulation domain.Model parameters are set per default values in Table 2.

Grid aggregated output
Grid aggregated results from both TopoSUB and distributed runs are analysed as deviations from BASE simulation for increasing sample numbers (TopoSUB: samples, distributed: pixels) that represent computational cost (Fig. 7).The simulated TVs approximate the BASE results well requiring 10 3 -10 4 times less computations.The convergence of results for the tested TVs at 100-200 samples in the sub-grid scheme suggests that we are able to reach a stable performance.A convergence of results with resolution is not observed in distributed simulations (except for the simple variable of air temperature) underscoring the importance of attention to scaling issues.Figure 7 also shows the improvement in aggregated information between a grid average computation (1 sample), as is common in climate models, and a 200 sample TopoSUB simulation, that is capable of approximating the BASE simulation (which is explicitly modelling sub-grid processes).

Statistical description of sub-grid behaviour
CDFs are calculated using aggregated sample weights and provide a rapid description of sub-grid behaviour by presenting the distribution of simulated TVs without the need for spatialisation.This technique enables rapid assessments such as percentage of permafrost in the simulation domain or the total SWE.A good fit is seen in all cases (Fig. 8).These statistics can be readily presented against topographic attributes to give more detailed understanding of the sub-grid e.g., permafrost extent by elevation band or exposures.

Spatialised output
Spatialised results are obtained by distributing LSM results based on the crisp or fuzzy membership of each pixel, resulting in a 25 m resolution mapping of TVs.Based on convergence of the NRMSE and correlation coefficient of the TopoSUB scheme it can be seen that, at least in our test case, the majority of performance is gained until 64 samples, after 258 samples a stable performance is achieved (Fig. 9).This represents a reduction of computational effort of three to four orders of magnitude compared with the BASE simulation (depending on required quality level), a similar result to Sect.5.1.Figure 10 gives density scatter plots of all TVs for 258 samples.TopoSUB is able to reproduce the BASE simulation with an NRMSE of 12-28 % depending on the TV.
Figure 11 provides a visual comparison of the simulation results for GST presented as deviation from BASE simulation (BASE TopoSUB), as well as a histogram of error distribution.The source of this error was investigated through a regression analysis of difference against PREDs.By restricting the dataset to values > 1 • C and < −1 • C we could ensure the signal was not masked by the (vast majority of) low error values.The model explained 47 % of variance (increased 1.0 samples correlation coefficient, r q q q q q q q q q GST SWE SWin Tair q q q q q q q q q q GST SWE SWin Tair Fig. 9.The predictive power of the lumped scheme is assessed using the (a) r-value (b) NRMSE at resolutions 2-400 samples with respect to baseline distributed grid simulation.The majority of performance is gained until 64 samples (first dotted line).After 258 samples a reasonably stable performance is achieved (second dotted line).
to 62 % by including interacting effects).A relative importance metric derived from decomposed r 2 value (Genizi, 1993) gives the percentage of variance explained by model attributable to each PRED, as follows: sin(aspect) = 38 %, cos(aspect) = 24 %, slope = 20 %, elevation = 14 % and sky view factor = 4 %.This shows that the spatial component of  the error is a reasonably even result of interactions among the PREDs, with only sky view factor being insignificant.

Temporal errors
Figure 12 shows the temporal error signature of the four tested target variables as 5-day mean values for the whole simulation domain, as absolute values for BASE, TopoSUB and difference between BASE simulation and TopoSUB.This figure shows that there is a temporal signature of the error of TopoSUB (at this temporal resolution), but also, that it is relatively small.Swin and SWE show the most obvious temporal trend with a stronger negative bias in SWin during winter months (approx. 2 W m −2 ) and an increasing positive bias in SWE during the main snow melt months (April-June), up to 7 mm.The temporal dimension is an extension of the multivariate problem and clustering may need to be tuned to fit certain seasons in much the same way as simulating different target variables (through informed sampling).

Model stability
Figure 13 shows results of 40 simulations of TopoSUB at each of four resolutions: 25, 50, 100, 200 samples.The purpose being to investigate stability of the tool and any resolution dependency of its repeatability.The deviation of each simulation from mean values of mean and quantiles 25/75 of all 40 simulations indicates reasonable repeatability even at low resolutions, as indicated by a low absolute deviation.This demonstrates that the result is not significantly sensitive to variations in the K-means clustering algorithm (which are generally small and diminishing with increasing sample number).

TopoSUB configurations
Figure 14 gives NRMSE for SUB (128 sample) simulations with the configurations: "fuzzy informed", "fuzzy simple", "crisp informed", "crisp simple".All TV results improve with fuzzy membership.All TV results except SWin improve with informed scaling.This is because all TVs are heavily influenced by PRED "elevation", except for SWin.During calculation of β mean SWin specific PREDs (e.g., sine of aspect) are down-weighted in favour of elevation -causing reduced performance for SWin results.This can be corrected for using the parameter W TV (Eq.5) depending upon user application.These results demonstrate that informed scaling and fuzzy membership are both useful tools to be applied in the pre-processor stage.Time-cost associated with these methods is incurred only once during pre-processing and benefit (improved performance) is received multiple times with LSM output.

Conclusions
This paper has described the TopoSUB scheme and has evaluated its performance by comparison with a high resolution distributed model (BASE), as well as given first indication of parameter choice.We have obtained promising results for the use of TopoSUB as robust and efficient tool to support the large area numerical simulation of processes in complex mountain topography.The informed sampling procedure is thought to be an appropriate method by which to sample a multidimensional sub-grid space without a priori knowledge of the relative importance of dimensions, with respect to simulated TVs.In the presented case, a sample number of approximately 64 samples proved sufficient to describe the wide range of elevations, slope expositions, slope angle and horizon elevations present in the study area, while 258 samples approximated well the results of a distributed 2-D simulation of 25 m resolution.While spatialisation proves a valuable tool for site specific studies, CDFs of sample results provide a rapid assessment of sub-grid behaviour.
There are obvious limitations which must be acknowledged due to the tool being based on a 1-D configuration, namely lateral mass transfers (such as snow redistribution, although this is commonly neglected even in distributed models) and flow modelling (surface/subsurface); processes which are not easily represented in a 1-D model.Additionally, at this stage we do not assume results derived from our specific test case (as defined by model, parameters and domain) will necessarily generalise, this would require more testing.
The ability of TopoSUB to efficiently approximate a distributed model as a reduced series of 1-D samples has several applications: 1. Multiple repeat simulations can be performed which facilitates exploring of uncertainty in initial and boundary conditions, model physics, parameterisations/parameter choice and also future climate scenarios.
2. Long temporal scales can be simulated allowing for modelling of transient effects.
3. Large areas can be simulated allowing for large scale assessment of current or possible future conditions.
www.geosci-model-dev.net/5/1245/2012/Geosci.Model Dev., 5, 1245-1257, 2012 4. TopoSUB is able to approximate the results of of a distributed grid that have been aggregated at coarse grid level.This suggests interesting prospects for the improvement of the representation of mountains or other environments with important fine-scale processes in coarse-grid models.
5. This approach is suitable for grid-computing infrastructure that is becoming increasingly common in many disciplines and, thus, allows flexible scaling to large tasks.
We envisage TopoSUB to be a useful tool in a wide range of numerical modelling applications in complex terrain, due to flexible choice of inputs, numerical models and output options.

Fig. 2 .
Fig. 2. The effect of maximum allowed iterations (iter.max)and number of random starts (nstart) of Kmeans on total within sum of squares (WSS) of resulting clusters (k = 128).
Fig. 4. (a) Visualisation of 128 samples generated by TopoSUB.(b) Polar plot showing the distribution of samples in terms of predictors used: elevation, aspect, slope and sky view factor.

Fig. 5 .
Fig. 5. Scheme of pre-processor output options, (A) no spatial, sample matrix of sample weights and environmental characteristics; (B1) spatial with crisp membership, single layer maps; and, (B2) spatial with fuzzy membership, multi-layer maps with a number of chosen fuzzy membership dimensions.

Fig. 7 .
Fig. 7. Aggregated statistics: mean (bold), 25th and 75th percentile (non-bold) of the sub-grid scheme at resolutions 1-1024 samples (blue) and distributed simulations (red) spanning resolutions of 4-10 6 pixels.Vertical lines indicate 16 and 128 samples.A stable performance is reached after the 128 sample level in all cases.Topo-SUB is able to approximate aggregated 2-D simulation at 10 4 less computations.Note logarithmic X-scale.

Fig. 8 .
Fig. 8. CDFs of mean annual simulation results derived from the sub-grid scheme (blue) based on 258 samples and a distributed simulation (red) based on 10 6 pixels.A good fit is reported by a KS-test (D).Aggregated summary statistics can be computed directly from the CDF, for example percentage area with MAGST < 0 • C (red dotted line, top-left figure).

Fig. 10 .
Fig. 10.Density scatter plot of TopoSUB/BASE after informed scaling and fuzzy spatialisation at 258 samples.Data presented is mean annual value for each pixel in the simulation domain.All TVs are reproduced with low error as reported by the correlation coefficient (r) and RMSE (computed over 10 6 pixels).The diagonal line represents y = x.

Fig. 11 .
Fig. 11.(a) A visual comparison of the simulation results for GST (128 samples) presented as deviation from BASE simulation (BASE TopoSUB) and, (b) a histogram of error distribution.Errors statistics show the error to be reasonable: RMSE = 0.6, bias = −0.15,standard deviation = 0.58.

Fig. 12 .
Fig. 12. Temporal characteristics of the mean domain values plotted at 5-day intervals for each target variable as simulated by BASE (red), TopoSUB (blue) and difference, BASE TopoSUB (green).Errors are relatively small in all cases.A stronger negative bias in SWin during winter months (approx. 2 W m −2 ) and an increasing positive bias in SWE during the main snow melt months (April-June) up to 7 mm.

Fig. 13 .
Fig. 13.Repeatibility of TopoSUB clustering success expressed as deviation from mean values of mean and quantiles 25/75 for 40 runs.Results at resolutions of 25-200 samples indicate reasonable stability even at low resolutions.A significant increase in stability is seen between 25-100 samples in all variables tested.

Table 1 .
Speeds t (min) of K-means in terms of parameters, cluster number k, pixel number p, and nstart n.