Spatio-Temporal Variations and Uncertainty in Land Surface Modelling for High Latitudes: Univariate Response Analysis

. A range of applications analysing the impact of environmental changes due to climate change, e.g. geographical spread of climate sensitive infections (CSIs), agriculture crop modelling, etc. , make use of Land Surface Modelling (LSM) to predict future land surface conditions. There are 5 multiple LSMs to choose from that account for land processes in different ways and this may introduce predictive uncertainty when LSM outputs are used as inputs to inform a given application. For useful predictions for a speciﬁc application, one must therefore understand the inherent uncertain-10 ties in the LSMs and the variations between them, as well as uncertainties arising from variation in the climate data driving the LSMs. This requires methods to analyse multivariate spatio-temporal variations and differences. A methodology is proposed based on multi-way data analysis, which extends 15 Singular Value Decomposition (SVD) to multi-dimensional tables, and provides spatio-temporal descriptions of agreements and disagreements between LSMs for both historical simulations and future predictions. The application underlying this paper is prediction of how climate change will affect 20 the spread of CSIs in the Fenno-Scandinavian and north-west Russian regions, and the approach is explored by comparing Net Primary Production (NPP) estimates over the period 1998-2013 from versions of leading LSMs (JULES,

ORCHIDEE-MICT respectively differ from the other LSMs. Differences between the climate forcing GCMs had a marginal effect up to 6% on NPP predictions out to 2100 without specific spatio-temporal GCM interaction.

Introduction
The rise in surface temperatures under global warming is predicted to be most severe in the Arctic, where it is already altering surface conditions and perturbing ecological systems (Overland et al., 2014). This will have multiple societal impacts, not 5 least on the health of animals and humans (IPCC AR5 WG2 A, 2014). The term climate sensitive infection (CSI) refers to diseases whose epidemiological aspects are driven, at least in part, by climatic factors (McMichael et al., 2006;Ebi et al., 2017;Cayol et al., 2017). In the Arctic, climate change is likely to cause enhanced CSI risk in terms of increased incidence, more frequent outbreaks, geographic spread of existing affected zones, and occurrence of newly affected zones (Pauchard et al., 2016;Sajanti et al., 2017;Waits et al., 2018) The complex ecology of CSI organisms presents a challenge to modelling and 10 predicting their epidemiology (Ostfeld, 2010;Carvalho et al., 2014;Ruscio et al., 2015;Sormunen et al., 2016;Li et al., 2016;White et al., 2018). However, such modelling is needed as disease vectors, such as ticks, mosquitoes, badgers and roedeer, which are associated, for example, with Lyme disease (Borreliosis) and Tularemia are expanding their ranges northwards (Jaenson et al., 2012;Jore et al., 2014;Andersen and Davis, 2017;Laaksonen et al., 2017;Blomgren et al., 2018).
Under climate change and human-induced land use changes, fragmentation of the landscape was found to affect Lyme disease 15 (Simon et al., 2014), whilst mosquito abundance was associated with outbreaks of Tularemia in boreal forest (Rydén et al., 2012). CSIs can also be non-vector diseases, as climate change may force increasing proximity and contacts between animals, e.g. pestivirus affects mammals (livestock or wild) and so reindeer (Kautto et al., 2012). This could threaten the successful bovine pestivirus eradication programs existing in Scandinavia since the 1990's (Tryland, 2013).
The Nordic Centre of Excellence (NCoE) CLINF, "Climate change effects on the epidemiology of infectious diseases and the 20 impacts on Northern societies" (www.clinf.org), is an interdisciplinary project supported by NordForsk (www.nordforsk.org), covering an area extending across Norway, Sweden, Finland and north-west Russia. Its aim is to study how climate change will affect the prevalence of human and animal CSIs and the consequences for Arctic societies. To do so it needs to characterise how a changing climate will change the environmental and societal conditions affecting a range of CSIs in Nordic regions. Besides predicting environmental changes likely to affect the spread of CSIs, CLINF also aims to gather and generate information on 25 the societal impacts of climate change. Achieving this aim requires tools to model land surface and aquatic changes under climate forcing. This paper focuses on land surface models (LSMs) and the extent to which existing LSMs could provide forecasts useful for the purposes of predicting CSI epidemiology.
An important factor in discussing the predictive value of these models is the variability in their outputs. This variability arises from two sources: variability in the climate drivers, since there are many Global Circulation Models (GCMs), and differences 30 between LSMs, whose core concepts are similar but with many differences in process representation and parameterisation.
This leads to three key questions: (i) How does the choice of the GCM affect the CSI-relevant outputs of a given LSM?
(ii) For a given GCM, how different are the CSI-relevant outputs of the different LSMs?
(iii) How do the joint effects of GCM and LSM differences translate into variability in predictions of CSI-relevant quantities?
Addressing these questions requires methods to describe spatio-temporal differences in models, and the first part of this paper describes such methods. The treatment here is relevant to a range of applications and is generic, but the evaluation of the methods in the latter part of the paper is couched in terms of differences between LSM predictions of Net Primary Productivity 5 (NPP), i.e. a single model output variable indicating vegetation activity, hence with relevance to CSI modelling involving changes in habitat for specific vectors, as well as carbon fluxes and ecosystem functioning (Koca et al., 2006;Rafique et al., 2016).
It is important that we quantify the uncertainty in any variable derived from an LSM model as a predictor in CSI modelling, so that the full uncertainty in the predictions (and associated risk) is available to public health decision-making. Typically, the 10 uncertainty in the predictions from a single LSM is poorly known, and we instead treat the spread in data simulated by a range of leading LSMs as a proxy for this uncertainty. Since Arctic CSIs are the underlying motivation for this work, we only consider LSMs that represent characteristics of Nordic areas, including high latitude processes, vegetation and landscapes. These are CLM5 (the Community Landscape Model version 5) (Lawrence et al., 2019); JULES (the Joint UK Land Environment Simulator) (Clark et al., 2011;Comyn-Platt et al., 2018); and two versions of ORCHIDEE (ORganizing Carbon and Hydrology in 15 Dynamic EcosystEms), ORCHIDEE-MICT (OR_MICT) (Guimberteau et al., 2018) and ORCHIDEE-HLveg (OR_HL) (Druel et al., 2017). The simulated climate data cover the historical period from December 1997 to December 2013, while for JULES we also analysed data from 100-year forecasts to the end of the 21 st century under forcing by 34 different GCMs (Comyn-Platt et al., 2018). The specifics of the four models and the driving climate data are briefly described in section 1.2.
Section 2 motivates the use of a multi-way methodology to characterise variations between LSMs, and the essentials of such 20 a methodology are described in Section 3. In Section 4 we use this methodology to analyse the differences between the four selected LSMs, while Section 5 shows how the methodology can be applied directly to differences between the LSMs. The same approach is then used in Section 6 to assess how the choice of a particular GCM affects the NPP predictions from the JULES LSM. Section 7 gives our discussion and conclusions.
1.1 LSM and ecological modelling aspects relevant to CSI prediction 25 Climate change is driving the spread of a range of CSI disease vectors (Zuliani et al., 2015;Andersen and Davis, 2017;Blomgren et al., 2018), so understanding the spatio-temporal distribution and evolution of characteristics, such as habitat suitability of these vectors or reservoirs, is essential. These characteristics can then be used in an ecological model that could be coupled with epidemiological models to estimate future risks of disease incidence, e.g. where and when the transmission risks are likely to be highest. For example, changes in abundance and extent of habitat suitability are important factors to be 30 considered in dynamic landscape epidemiology modelling (Lambin et al., 2010). Changes due to global warming could affect both abundance and geographical extent, and also extend the period of transmission, e.g. by affecting the vector life cycle (Rose et al., 2015). Extreme weather conditions and events may either introduce outbreaks of abundance, thereby increasing the risk of pathogen occurrence in disease vectors, or wipe out a species at a given location.
Under different scenarios of climate drivers, such as the Representative Concentration Pathways (RCPs) developed by the Intergovernmental Panel on Climate Change (IPCC), LSMs can simulate future atmospheric and land conditions that can be related to vector habitat suitability. Variables simulated by combined climate and land surface models, such as surface 5 temperature, soil moisture, precipitation and land cover, can be used in ecological models or as part of an epidemiological model, e.g. for a species distribution model (Booth et al., 2014). However, the predictive uncertainty in these variables may lead to significant uncertainty in the predictions from CSI modelling (Asghar et al., 2016). This paper describes and quantifies the spatio-temporal uncertainty arising from the choice of LSM alone, i.e. without assessing its impact on CSI predictions, but provides an essential component in understanding the uncertainty in any statistical or mathematical predictions of CSI 10 epidemiology and ecology (Beale and Lennon, 2012;Zuliani et al., 2015;Metcalf et al., 2017) that use LSM outputs as predicting variables.

Land Surface Model and Data Description
The four LSMs used in this study, CLM5, two versions of ORCHIDEE (OR_MICT and OR_HL) and JULES, were chosen because of their high degree of maturity and their ability to model characteristics of Nordic areas, including high latitude pro-15 cesses, vegetation and landscapes. Table 1 summarises these characteristics; details can be found in the references. OR_MICT (Guimberteau et al., 2018) includes major high latitude adaptations, including snow and soil thermal interaction, plant primary productivity constrained by high latitude conditions, and soil carbon stocks with feedback dynamics. OR_HL (Druel et al., 2017(Druel et al., , 2019 adapts ORCHIDEE with specific plant functional types (PFTs) such as non-vascular plants (mosses, lichens), Arctic C3 grass and boreal shrubs. CLM5 (Lawrence et al., 2019) includes permafrost modelling and snow processes, C3 20 Arctic grass and deciduous boreal shrubs as part of its 15 PFTs (see Appendix B) but no non-vascular plants. The version of JULES (Clark et al., 2011) used here has been extended to be suitable for high latitudes (Comyn-Platt et al., 2018) by including processes such as permafrost-carbon feedbacks (Burke et al., 2017).
For all the LSMs, the initial PFTs were derived from land cover maps. JULES and the two versions of ORCHIDEE use the land cover product from the European Space Agency Climate Change Initiative (ESA-CCI) (Poulter et al., 2015). The 25 supplementary material to Druel et al. (2017) describes the correspondence between land cover and the added Arctic PFTs.
CLM5 uses the Land Use Harmonised data version 2, a product of the Land Use Intercomparison Project (LUMIP) (Lawrence et al., 2016) to define its initial spatial distribution of PFTs. For the historical analyses, the data were re-gridded to the finest grid-spacing, 0.5°E × 0.5°N, by simple disaggregation which introduces a limitation when comparing the LSMs. All analyses were performed for a sub-area of the CLINF zone between 4.5°E-34.5°E and 58.5°N-70.5°N. Note that the climate forcing 30 data are not the same for the different LSMs (see Table 1) since the LSM data were provided by different modelling groups, each of which uses preferred GCMs. This is unlikely to have any significant impact on the LSM comparisons (see Section 5). the LSMs provide NPP for each PFT, the PFT dimension could also be added, but this is not done here.
Analysing such structured datasets to understand spatial, temporal and between-model variations can be challenging when there are long-tail distributions (as is the case in our dataset: see Fig.1, which shows the histogram of NPP values in the combined historical datasets simulated by the 4 LSMs) which preclude the use of classical geostatistical methods, and due to their multi-variate :::::::::: multivariate nature. For 2-way tables, Singular Value Decomposition (SVD) is a powerful tool to extract associations of variables and patterns within data, e.g. clusters and trends. The SVD of a data matrix finds pairs of vectors (components) that successively extract decreasing fractions of the variation in the data and are uncorrelated with previous pairs. Visual description of these optimal vectors can be obtained by plotting the component weights, e.g. a Spatial effect as a 5 map and a T emporal effect as a time series plot.

From Singular Value Decomposition to multi-way data analysis
Let X be a n × p matrix, which we can regard as a collection of n p-dimensional vectors or p n-dimensional vectors. The matrix X t X is positive semi-definite, so all its eigenvalues are positive, and its eigenvectors, ϕ h , are mutually orthogonal, i.e.
The matrices X t X and XX t have the same eigenvalues, σ 2 h , and the sum of squares of the elements of 20 X is given by x is the matrix X vectorised as a np-dimensional vector. The SVD of any matrix X is defined by the series of decreasing σ h , the singular values, each associated with a pair of unit vectors ϕ h (p-dimensional) and ψ h (n-dimensional), with ψ t h ψ h = ϕ t h ϕ h = 1, which explain a fraction σ 2 h /( h σ 2 h ) of the variability of X (defined as the sum of squares of the elements of X), i.e. ϕ t h X t Xϕ h = σ 2 h = ψ t h XX t ψ h . Hence SVD can be used for dimension reduction by defining a p -dimensional subspace (p < p) that captures most of the variability in X: For a suitable p , the residual variation = h>p σ h ψ h ϕ t h is small enough to be considered as insignificant. As shown in equation (1), this decomposition can be written in tensorial form, since ψ h ϕ t h = ψ h ⊗ ϕ h . The rank-1 matrix ψ h ϕ t h is known as a decomposed rank-1 tensor (Leibovici, 2010). The term σ 1 ψ 1 ⊗ ϕ 1 is the best rank-1 tensor approximation to ::: the :::::: matrix X in the sense of capturing the maximum fraction of variability in X among all rank-1 tensors, i.e. the maximum value of σ = ψ t Xϕ. Subsequent rank-1 tensors in the decomposition in equation (1), given by the other eigenvectors, are orthogonal to the previous ones and successively extract decreasing fractions of the variability. Matrices can be seen as order 2 tensors and 10 multi-way tables as order k tensors, where k is the number of dimensions of the table. For k = 2 the SVD can be seen as an optimal basis vector system in each dimension and in the tensor space, and generalisations of SVD to tensors of order k ≥ 3 aim to find equivalent optimal systems. Several algorithms to extend SVD to tables with more than 2 entries have been proposed (Tucker, 1966;Carroll and Chang, 1970;Harshman, 1970;Kroonenberg, 1983;Leibovici et al., 2007;Leibovici, 2010) and development of methods and their 15 applications is is still very active (Demšar et al., 2013;Kroonenberg, 2016;Takeuchi et al., 2017;Lock and Li, 2018). Most extensions aim to find an optimum decomposition of a multi-way table that allows dimension reduction by looking for a decomposition similar to equation (1) under specific optimisation criteria. ::: The :::::::: algebraic ::::::::: embedding ::: of ::::::: two-way :::: data :::::: tables :: as :::::: tensors ::: of :::: order :: 2 :::::::::: (equivalent :: to :::::::: matrices) :::: and ::: by ::::::::: extension, :: of :::::: k-way :::: data ::::: tables ::: as :::::: tensors ::: of ::::: order :: k, ::::::::: facilitates ::: the ::::::::::: understanding :: of ::::: these :::::::::: extensions. For a multi-way table X with k ≥ 3 entries this takes the generic form: where the h i index the basis vectors of the individual vector spaces making up the k-dimensional data table and expresses the residual of the approximation given by the summation. This residual depends on the method and the number of components used in the decomposition, and can be zero (as would be the case if we retained all the terms in a SVD). 25 The decompositions carried out by the CANDECOMP and PARAFAC methods (Carroll and Chang, 1970;Harshman, 1970) fix the number of rank-1 tensors in the decomposition but do not impose an orthogonality constraint, while PCA-3modes (Kroonenberg, 1983) considers both orthogonality and rank within each vector space. However, all three methods need to choose in advance the number of rank-1 tensors in their optimisation and obtain decompositions that are not nested as with SVD, in which the rank p approximation of X contains the approximation obtained for p (with p > p ). This property is 30 often desirable for environmental data analysis (Frelat et al., 2017), as decomposition of the variance or sum of squares has a practical interpretation.
To address this, Leibovici and Sabatier (1998) developed the PTAk method, which is a hierarchical decomposition giving nested approximation by construction. For k = 2, the PTAk algorithm is the same as SVD, while for k = 3 it is given by: The notation ⊗ i means that the vector on the left takes the ith place in the tensor product, e.g. ϕ 1 ⊗ 2 (α ⊗ β) = α ⊗ ϕ 1 ⊗ β and 10 ".." indicates the contraction operation (defined in Appendix A along with definitions of the other notation used in equation (3)). Note that the PTA3 algorithm is recursive as the last line of equation (3) calls another PTA3. This process can be continued until it leads to a null table, but normally a stopping rule is imposed by requiring the decomposition to capture a prescribed fraction of the overall variability or specifying the desired number of order k PTs ::::: rank-1 ::::::: tensors ::::::::: (right-hand :::: side :: of :::::::: equation :: 3).

15
Similarly to SVD, the PTA3 algorithm first retrieves the rank-1 tensor approximation to X, σ 1 ψ 1 ⊗ ϕ 1 ⊗ φ 1 , that captures the maximum possible fraction of the variability in X ::: and :: is :::::: termed ::: the :::: first ::::::: principal :::::: tensor :::: (PT) :: of :: X. The second, third and fourth lines in equation (3) correspond to optimisations associated with this first PT in which the decomposed tensors share one of the components in the the first PT. The corresponding PTA2 analyses are complete SVD decompositions into series of tensor products. Given this decomposition, descriptive statistics or plots of the triple of components (ψ 1 , ϕ 1 , φ 1 ) can then be 20 used to visualise the pattern or effect associated with the fraction of the variability captured by each of the tensors.
The generalisation of equation (3) to k-way data tables is straightforward. In a PTAk decomposition, the first rank-1 tensor will have associated PTA(k − 1)'s which will recursively end up at associated PTA2's, i.e. SVDs.

Spatio-temporal variations of NPP across the LSMs
This principal aims of this section are to perform a PTA3 analysis of the 3-way Spatial × T emporal × LSM table X of NPP 25 and to interpret the results. However, it is useful to first examine some of the properties of the distributions of NPP for each LSM. The histogram of the NPP values in the full data table X, displayed in Fig.1, conceals distinct differences between the LSMs. Some of these differences are indicated by Table 2, which gives the mean NPP and sum of the squares of NPP for each LSM, and Table 3, which shows for each LSM the fraction of NPP values in each decile of the reference distribution in Fig.1.
In Table 2, OR_MICT stands out by its low mean NPP (23% less than JULES) and low variability (significantly less than the 30 other LSMs, and 37% less than OR_HL). The LSMs also exhibit different distributions (see Table 3): notably CLM5 has 35% of its NPP values in the first decile of the reference distribution, while OR_HL and JULES have very few values in this decile, and the decile with peak occupancy is different for all four LSMs. However, all the LSMs place around 10% of their NPP estimates in each of the higher deciles (70% to 100%). These distributional differences tell us nothing about the spatio-temporal differences between the LSMs, and for that we use the decomposition provided by the R package PTAk (Leibovici, 2010) of which the first ten terms are displayed in Fig.2.
This describes the hierarchical and nested decomposition of the sum of squares of X into PTs and associated PTs. Each row 5 corresponds to a PT, identified by a label and number, -no-, and its singular value, e.g. vs111 and -no-1 correspond to the first line of equation (3) giving the best rank-1 approximation of X, with singular value σ 1 = 2.7147e − 05. The row with label vs222 gives the singular value corresponding to the next order-3 rank-1 approximation, which corresponds to the recursive step in the last line of equation (3).
The rows between vs111 and vs222 correspond to PTs associated with vs111, which are derived from PTA2s, i.e. SVDs.

10
The labels given to these decomposed components start with the dimension of the component used in contracting the tensor X (see Appendix A) and continue with the label of the PT from which they are derived and the dimensions of the contracted tensor, e.g. 1152 vs111 193 4 identifies the results from the PTA2 of the 193 × 4 matrix X..ψ 1 (i.e. an SVD), where ψ 1 is the 1152-dimensional vector forming the Spatial component of PT -no-1. Therefore the associated PTs -no-3 and 4 have the same Spatial component as tensor -no-1. Similarly, the rank-1 tensors -no-6 and 7 are associated PTs along the Table 3. Deciles (q) of the reference NPP distribution given by Fig.1 and the percentage of NPP values in each decile observed for each LSM.
An LSM whose NPP values had a distribution similar to the reference would have 10% in each decile; entries in bold indicate departures of more than 2% from 10%.
q -1.25e-09 -5.15e-11 4.58e-11 1.30e-09 5.67e-09 1.53e-08 2.65e-08 4.14e-08 5.80e-08 1.11e-07 OR_MICT  The contribution by the main PTs decreases from vs111, vs222, vs333, etc. Each of the associated tensors makes a smaller contribution than its main PT but this may be larger than the next main PT, e.g. tensor -no-3 captures more variability than tensor -no-11. There is no particular ordering in the tensors associated with different components, between -no-3 which is associated with the Spatial component and -no-6 which is associated with the T emporal component, but the PTs associated with a given component are ordered since they derive from the same PTA2 (i.e. SVD), e.g. -no-3 precedes -no-4. Fig.2 then allows one to select the PTs or associated PTs that successively capture the variability in X.
It is helpful to visualise the first PT, whose components are displayed in Fig.3, as an optimal approximation to the initial 1152 × 193 × 4 data table in which each of the 4 layers is the same 2-D spatio-temporal "map", but scaled by the weight for a particular LSM, given by φ 1 . The spatial pattern at each time is the same (ψ 1 , as in Fig.3(a)), but with a weight appropriate 5 to that particular time. Similarly, the time series at each spatial location is the same (ϕ 1 , shown as Fig.3(b)), but with a weight appropriate to that location. To recover the NPP from this approximation at a particular position, time and for a given LSM, the corresponding values in ψ 1 , ϕ 1 and φ 1 are multiplied together and then multiplied by its singulat ::::::: singular value, σ 1 . Exactly the same construction applies to each of the rank-1 tensors in the decomposition.
The Spatial effect (Fig.3(a)), which is always positive, places higher weights in Sweden, the Baltic states and north-west Since these spatial and temporal patterns are the same for all the LSMs, the difference between them is expressed by the LSM weights (Fig.3(c)). For identical LSM simulations, the weights would be 1/2, since each vector in the decomposition is normalised to unity (i.e. φ 2 11 + φ 2 12 + φ 2 13 + φ 2 14 ) = 4φ 2 11 = 1), but the weights lie between 0.460 and 0.523, with JULES 25 and OR_HL respectively giving values 3% and 5% higher than for equal weights, and OR_MICT giving a value 8% lower.
Hence there is only a weak dependence on the LSM in this first PT. The proportion of the variability in the first PT due to each LSM is given by the squares of the LSM weights, i.e. 21.2%, 27.4%, 25.1% and 26.3% for OR_MICT, OR_HL, CLM5 and JULES, respectively. Multiplying these values by σ 2 1 gives the sum of squares of NPP in the spatio-temporal maps for vs111 for each LSM (see Table.2). Several points should be noted about the approximation to X given by vs111:  Table 2) are in the ratio 1 : 1.37 : 1.25 : 1.29. Hence the first PT correctly picks up the ordering of the variability amongst the LSMs, but not its full value, since it is effectively a smoothing of the dataset.  3. The mean NPP represented by vs111 is 1.834 × 10 −8 kg m −2 s −1 , which is very close to that of the mean of X (1.824 × 10 −8 kg m −2 s −1 ), though the individual NPP spatio-temporal maps for each LSM track the original mean 5 NPP less closely (+3.6%, -0.7%, +6.3% and -5.6% for OR_MICT, OR_HL, CLM5 and JULES respectively; see Table   2).
As noted above, recovering the NPP at a particular position, time and for a given LSM in vs111 requires multiplying together the corresponding weights in the Spatial, T emporal and LSM dimensions and then multiplying by the singular value. So, for example, the maximum value of NPP in the first PT over the whole time-period is in July 2013, in the darkest red cell of  The second best PT in the decomposition, -no-6, is a T emporal-associated PT, so has the same T emporal component as vs111, and expresses 3.72% of the variability. Its Spatial component ( Fig.4(a)) has positive (red) weights in the north and west and negative (green) weights to the south and east. The most striking feature of this tensor is in its LSM component ( Fig.4(c)) which shows a marked contrast between OR_HL, with a large negative weight, and the other LSMs, for which the weights are significantly smaller and positive. Hence, after multiplying the weights in the different components, all the LSMs 15 except OR_HL will see an increase in NPP in the red areas in the summer months and reduce it in the green areas, while the opposite effect occurs for OR_HL. When the T emporal weights are negative, as occurs for most of the winter, these sign changes in NPP are reversed. As can be seen from  Here the T emporal effect is positive for the months from August to October, close to zero for November and July, and negative for the other months, especially April to June. CLM5 and JULES have large positive and negative weights respectively while OR_MICT has a smaller negative weight and OR_HL has a weight which very close to zero. Hence for CLM5 this tensor acts to increase NPP from August to October and reduce it for all other months except November and July, while for JULES and OR_MICT it does the opposite. As expected from the weights, including this 25 tensor mainly acts to improve the fit of CLM5 and JULES to their original values (Table 2).
Principal tensor -no-9 is the fourth best in the decomposition and captures 1.38% of total variability. Since it is associated with the LSM component of PT -no-1 it is the same for all LSMs. Its Spatial component Fig.6(a)  in April to June it does the opposite. These effects will be slightly greater for OR_HL because of its greater weight. Though its contribution to the overall sum of squares is only 1.38%, it provides improvements for all LSMs (see Table 2).
None of the other PTs contributes more than 1% to the overall variability and their components are not displayed, although the contributions for all terms in Fig.2 are given in Table 2. For example, the next best PT (-no-11), which derives from the   Table 2).
Overall, this analysis shows that a single optimal spatio-temporal pattern, with slightly different weights for the four LSMs (up to 14% maximum difference), provides a reasonably good approximation to all their estimates of NPP, capturing between (a) Spatial dimension (b) Temporal dimension (c) LSM dimension Figure 5. Plots for PT -no-3 associated with PT -no-1 along the Spatial dimension, which is therefore identical to Fig.3(a)); it represents 1.72% of the variability.
87% and 93% of the variability in the individual models, as well as around 90% of the variability in the combined LSM dataset.
The next best adjustment to this pattern is a spatial correction that principally applies to OR_HL and significantly improves the fit of the approximation to this LSM, with only small improvements for the other LSMs. The second best adjustment adds a temporal pattern that mainly affects CLM5 and JULES and improves the fit to these LSMs, with less effect on OR_MICT (a) Spatial dimension (b) Temporal dimension (c) LSM dimension Figure 6. Plots for PT -no-9 associated with PT -no-1 along the LSM dimension, which is therefore identical to Fig.3(c)); it represents 1.38% of the variability.
and none on OR_HL. The third best adjustment adds a new spatio-temporal pattern whose spatial component is roughly the opposite of that in the first PT (i.e. it is spatially similar but with opposite signs) but a quite different temporal component that is positive in the later summer months, negative in the late spring and early summer months, and roughly zero at other times.
. The improvement in the overall fit from the next best PT and all succeeding ones is less than 0.9%, and, although in two instances the fits to individual models improve by over 1%, in most cases the improvements are much smaller (see Table 2).
Summing the 10 PTs whose individual contribution to the overall variability exceeds 0.1% (Fig.2) provides an approximation to the overall data table that captures 98.4% of the overall variability and between 97.4% and 99.0% of the variability in the individual LSMs (Table 2). However, also of interest is the point-wise goodness of fit of the approximation, not just the variability it captures. This is represented by the table of residuals, i.e. the term in eq. (2). Around 75% of the absolute values of these residuals are less than 8.4% of the overall mean NPP, so in most cases there is a good point-wise fit to the original 5 data, but the maximum absolute value of the residuals (4.83 × 10 −8 ) is around the third quartile of NPP (3.44 × 10 −8 ). Hence, in some cases the approximation may be significantly different from the correct value despite the residuals contributing less than 1.62% to the overall variability.

Analysing differences between the LSMs
Section 4 identified differences between the LSMs captured by an optimal decomposition of the associated 3-way table. In this 10 section we instead directly analyse the variability in the differences between the LSMs, in order to localise where and when the LSMs disagree and thus to quantify spatio-temporally the uncertainty in NPP associated with the choice of a particular LSM.
We   Fig.7 displays the histograms of (a) the absolute values of normalised differences, which has a peak near zero but also a fairly long right hand tail, and (b) the absolute values of (NPP 1 -NPP 2 )/(NPP 1 + NPP 2 ) which is fairly flat across most of the range, with a small peak near zero, but with a large peak peak near 1. The latter indicates that for many times and places the 20 NPP values in one LSM are very small relative to one of the other LSMs. This occurs much more frequently in winter when CLM5 gives NPP values that are very small compared to those from the other LSMs. However, since NPP is small in winter, these large relative differences have little impact on overall annual production. Indeed Table 2 shows that the mean annual NPP from CLM5 exceeds that from OR_MICT.
The results of the PTA3 for the 1152 × 193 × 6 table of normalised NPP differences are shown in Fig.8. The first and second PTs respectively extract 43.4% and 21.7% of the variation, both with well-structured patterns in their components. The first, shown in Fig.9, has a spatial pattern with negative (green) values areas to the south and east and positive (red) values in the 5 north and west, as well as south-east Finland. The T emporal component is always positive and displays a seasonal effect ( Fig.9(b)) with the same ordering of the months as Fig.3. All the differences involving OR_HL have significant weights but for the other differences they are close to zero. Hence the effects of this principal vector essentially translate into differences between OR_HL and the other LSMs. Taking into account the signs of the Spatial, T emporal and LSM weights (the last to be interpreted as an LSM difference), this means that for this PT over the whole time period CLM5 > JULES > OR_MICT > 10 OR_HL in the red areas in Fig.9(a) while these orderings are reversed in the green areas. However, the small weights on the differences not involving OR_HL indicate that the other LSMs all give similar NPP values. to differences in data resolution before grid transformation but also occurs where C3 grass is the dominant PFT (all LSMs).
All LSM differences have positive weights except CLM5 -JULES, which is negative but small, and all differences involving OR_MICT have significantly larger weights than other combinations. Since the temporal component is everywhere positive (Fig.9(b)), the net effect is that OR_MICT > OR_HL > JULES > CLM5 in the red areas (with a high value north of Lake Ladoga), and this order is reversed in the green areas. However, the differences between LSMs other than those involving OR_MICT are small.
The Spatial component of PT -no-7 (Fig.11(c)) is weakly positive except for a very small area near Tromsö in northern Norway. All the LSM differences involving JULES have negative weights and have greater magnitude than the other differ-ences, which are all positive, meaning that JULES > OR_HL > OR_MICT > CLM5 everywhere except near Tromsö, where this ordering is reversed.
As indicated by Fig.8, the next best PTs (not displayed) are -no-13, associated with the Spatial component of vs222, -no-9 and 10, associated with the LSM component of vs111, and -no-19, associated with the LSM component of vs222.
Hence PT -no-13 modulates the temporal pattern of differences depicted in Fig.10 with a distinct temporal pattern that has 5 different positive weights for each of the LSM differences (the contribution from OR_HL -CLM5 is almost zero and OR_MI -JULES gets the larger positive weight. A contrast between July (positive weights) and May (negative weights) stands out clearly from the other months by the size of its contribution to the variability, for reasons which are not clear. In Fig.10, July and and OR_MI -JULES weights were close to zero. Because PTs -no-9 and 10 are associated with the LSM component of vs111, the spatio-temporal table given by summing the Spatial x T emporal terms in all three PTs can be analysed together; this 10 would mainly reveal spatio-temporal differences between OR_HL and the other LSMs (see Fig.9(c). However, this combined analysis cannot be displayed as separate Spatial and T emporal plots. With the same LSM weights as in Fig.10, PT-no-19 exhibits a clear north-south gradient and a temporal pattern in which June clearly contributes more to the variability than the other months. This is similar to what is seen for July in PT -no-13, again for unknown reasons. All the rest of the PTs cumulatively contribute only 10% to the overall variability and individually less than 0.8%.

15
Also analysed was the variability in the quantity |NPP 1 −NPP 2 | / | NPP 1 + NPP 2 | but this is not displayed, since its main contribution is to show that the large peak near 1 seen in Fig.7 (plot (b)) can mainly be attributed to small values of CLM5 relative to the other LSMs in winter in the north of the region.
The analysis in this Section adds significantly to that in Section 4 by providing specific information on the times and places where the LSMs differ and by how much. However, in this case no single spatio-temporal pattern strongly dominates the 20 variability so interpretation of the analysis requires consideration of several such patterns. Nonetheless, the three best PTs capture around 76% of the variability in the LSM differences. The first essentially tells us that over a well-defined spatial pattern and a clearly ordered temporal pattern that with a maximum in summer and a minimum in winter, OR_HL gives different values from the other LSMs, which are all similar. The second PT principally identifies times and places where CLM5 differs from the other LSMs, while the third does the same for OR_MICT.

Climate forcing uncertainty
This section analyses the effects of different GCM drivers on the NPP estimated by JULES, so is a partial answer to question (i) in Section 1. Two global warming scenarios that stabilise at 1.5°C and 2.0°C above pre-industrial levels by year 2100 were used, with 34 GCMs as climate forcing in JULES (Comyn-Platt et al., 2018). The ensemble of the GCMs is taken to represent the uncertainty in climate prediction, from which one can get an idea of the associated uncertainty in the JULES estimates 30 of NPP. Note however, that this commonly-used approach to quantifying climate uncertainty is not entirely satisfactory, since it identifies inter-GCM model variability with the internal uncertainty in climate prediction (Hawkins and Sutton, 2009;Kay et al., 2015).
For each scenario a PTA3 analysis was performed on a Spatial × T emporal × GCM GCMs.
Over the 100 years, all months exhibit an initial increase, which is sharper for the 2°C scenario, followed by a flattening out and minor decrease; this decrease sets in around 2070 for the 1.5°C scenario and slightly later for the 2.0°C scenario. The maximum increase from 2000 (indicated on each monthly curve in Figs. 12(b) and 13(b)) is higher in every month for the 2.0°C scenario, e.g. 20% and 32% in July for the 1.5°C and 2.0°C case, respectively. The differences between the GCMs are indicated 10 by histograms of the relative deviation of the GCM weights from uniform weighting (Figs. 12(c) and 13(c)). These differences are up to 7% for the 2.0°C scenario and 4.5% for the 1.5°C scenario. For both scenarios, the groups of GCMs giving lowest or highest difference from equal weighting was the same, though the precise ordering was different (see Appendix C). If the singular value associated with this first PT are expressing the same amount of variability, the latter is 10% higher for the 2.0°C case than for 1.5°C, which simply expresses the sharper increase of NPP values produced under a more intense warming.
Global statistical differences were found between the LSMs, with OR_MICT exhibiting significantly lower mean NPP and variability than the other LSMs, and CLM5 producing a very high proportion of low values associated with the winter season, particularly in the north of the CLINF region. However, all the LSMs tend to agree for higher NPP values (above the 70% decile), which mainly indicates that they give similar values in summer. Despite these global differences, to a first 30 approximation the spatio-temporal behaviour of all the LSMs could be well-fitted by the tensor product of a single spatial and temporal pattern, in which the west and north of the region exhibited lower NPP values than the east and south, and there was a strong seasonal pattern. Differences between LSMs for this single pattern were fairly small, with weights lying between 92% and 105% of an uniform weighting of 0.5 or 14% maximum difference between them. This combined pattern captured around 90% of the overall variability in simulations covering 16 years for the whole Fenno-Scandinavian region. Across this time-period, this first approximation displayed statistically significant increases in NPP from May to September in , with the largest increase in the earlier months. This is likely to be caused by the growing season starting earlier and lasting longer.
The LSM requiring most adjustment to this first approximation for an improved fit was OR_HL; this adjustment is in the 5 spatial pattern, decreasing the spatial weights in Norway and northern Finland and increasing them in Sweden and southern Finland. The next adjustment, which has no effect on OR_HL, is to modify the temporal pattern; this particularly improves the fit to CLM5. The approximation achieved with just these adjustments captures 95% of the overall variation and between 93% and 96% of the variation in the individual models. It also has the advantage of being fairly simple to interpret because OR_HL dominates the first adjustment while CLM5 (and to a lesser extent JULES) dominates the second. Further terms in the 10 approximation yield smaller gains that tend to be spread more evenly across the LSMs.
While the first analysis provides information on temporal and spatial patterns characterising the LSMs, more specific information on how they differ is gained by analysing their differences. Here no single pattern dominates the overall variability between the LSMs, but the three best PTs capture around 76% of the variability in the LSM differences, and can be fairly well interpreted in terms of how individual LSMs differ in space and time from the others. Successively they show where and when 15 individually OR_HL, CLM5 and OR_MICT differ from the other LSMs, and also where different LSMs agree.
Our analysis of the impact of the choice of GCM on the simulations of NPP was restricted to runs with JULES out to 2100 driven by 34 different GCMs. This showed that a single spatio-temporal pattern captured over 99% of the variability of NPP in the combined dataset for climate change scenarios leading to either 1.5°C or 2.0°C atmospheric warming, and that none of the GCM weightings differed by more than 3% from uniform weighting (maximum difference of 6%). The temporal pattern 20 showed increases of NPP up to the 2070s, with small decreases thereafter. Although this analysis was only carried out for JULES, there is no reason to expect different findings for the other LSMs ::: (as ::: this :: is ::::: strong :::::: effect, :::: 99% :: of ::::::::: variability, :::::::: coherent :::: with :: the :::: first ::: PT :: in ::: the ::::::: analysis :::: with ::: the : 4 :::::: LSMs).
Returning to the three key questions posed in Section 1: (i) How does the choice of the GCM affect the CSI-relevant outputs of a given LSM? 25 (ii) For a given GCM, how different are the CSI-relevant outputs of the different LSMs?
(iii) How do the joint effects of GCM and LSM differences translate into variability in predictions of CSI-relevant quantities? the analysis in this paper suggests that, at least for NPP, we can neglect the effect of different GCMs and need only deal with question (ii). Quantitative answers are provided to this question both in terms of spatio-temporal patterns and differences and similarities of LSMs. However, we have only considered one of the six variables listed at the start of Section 2 that are 30 considered to be of major importance for Climate Sensitive Infections (CSIs), and may find different behaviour for the others.
In particular, initial investigation indicates very different representations of land cover between the four LSMs and how land cover will evolve under climate change in the 21st century. This variable is likely to be the one showing most differences between the LSMs because it is very much controlled by the PFTs used, how they are parametrised, and the rules by which PFTs compete over time.
Of significant interest would be analysis of multiple variables and their co-variation. We intended to address this issue in a future paper using the PTAk method used here, since this can be readily extended to multiple variables. While this does not present any methodological difficulties, it will only become clear how useful this is when we find how easy it is to interpret the 5 outputs of the analysis.
The next major step is to couple the findings from this paper (and its extension to other variables) to ecological models for CSI vectors and statistical epidemiological models in order to establish the sensitivity of predicted CSI behaviour under climate change to the choice of GCM and LSM. Currently only a small number of CSIs have well-developed predictive models (notably tularemia (Rydén et al., 2012;Andersen and Davis, 2017;Desvars-Larrive et al., 2017); Lyme disease (Simon et al.,10 2014; Li et al., 2016)) and these will provide the basis for such a study. However, CLINF is in process of developing more comprehensive statistical CSI models at high latitudes, which will lend themselves readily to combination with the approach adopted in this paper.
Code and data availability. The analyses were performed using the R package PTAk (https://CRAN.R-project.org/package=PTAk). The multi-way data tables used in the paper can be requested from the first author. CLM5.0 is publicly available through the Community Ter- Appendix A: Contraction operator and orthogonal projector

A1 Contraction
For X and Y two multi-way data tables n×p×q, their inner product is defined as < X, Y >= ijk X ijk Y ijk . The contraction 20 operation .. is the extension to tensors of the linear combination of the columns or rows of a matrix to give a vector. If X is a tensor of order 3, equivalent to a table n × p × q, then with the variables (u, v, w), vectors of length n, p and q, respectively, the contraction X..u is a p × q matrix with (X..u) jk = i X ijk u i , the contraction X..v is a n × q matrix with (X..v) ik = j X ijk v j , and X..w is a n × p matrix with (X..w) ij = k X ijk w k . Contacting X successively by two vectors gives for example (X..u)..v = ij X ijk u i v j = ij X ijk (u ⊗ v) ij = X..(u ⊗ v) and X..(u ⊗ v ⊗ w) is equivalent to the inner 25 product for the multi-way data tables.

A2 Orthogonal projector
Without loss of generality let u, v and w be unit vectors of dimensions n, p and q respectively. If X is a tensor represented by an n × p × q array, one can write X = (a ⊗ b ⊗ c)β + = P (a⊗b⊗c) X + P (a⊗b⊗c)⊥ X, where P a⊗b⊗c = (a ⊗ b ⊗ c)β is the linear orthogonal projection of X onto a ⊗ b ⊗ c and P (a⊗b⊗c)⊥ X = X − (a ⊗ b ⊗ c)β. From the orthogonality constrains, β = X..(a ⊗ b ⊗ c), so P a⊗b⊗c X = (a ⊗ b ⊗ c)X..(a ⊗ b ⊗ c).
Moreover, if X = (x ⊗ y ⊗ z) then P a⊗b⊗c X = P a x ⊗ P b y ⊗ P c z. This property extends easily to any subspace of E, F , and G, i.e P E1 ⊗ P F1 ⊗ P G1 is equivalent to P E1⊗F1⊗G1 . Appendix C: GCM's Weightings from the analysis in section 6

Appendix
The acronyms of the 34 GCM are derived from the information given in "  Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This work was carried out under NordForsk funding to CLINF, a Nordic Centre of Excellence (NCoE) led by Pro-