Interactive comment on “ Estimation of land-use change using a Bayesian data assimilation approach

Figure 4 should be improved. In the current figure the CI dominate the signal making it uninformative to show the prior and observations. If that is the message of this figure, then search for a more elegant way to show it (could be a table). The current presentation already uses different ranges of the Y-axis but even then for some rows of subplots the range is not completely used. We tried several ways of presenting these data, and we’re not sure there is a better


Introduction
Human-induced land-use change has a substantial impact on biodiversity and both biogeochemical and hydrological cycles (Post & Kwon, 2000;Gitz & Ciais, 2003;Levy et al., 2004;Newbold et al., 2015;Piano et al., 2017).The importance of representing it in models of the climate, hydrology, and ecosystem processes is increasingly recognised (Martin et al., 2017;Prestele et al., 2017;Quesada et al., 2017).However, although changes in land use tend to occur incrementally over small areas, data on land-use change are typically limited in spatial and temporal resolution (Alexander et al., 2017).Furthermore, changes in land use may be rotational or involve transitions between multiple land-use classes over time, such that the gross area undergoing land-use change may be much larger than the net change in area (Fuchs et al., 2015;Tomlinson et al., 2017).From the point of view of modelling ecosystem processes, it is these fine-scale gross changes that we need to represent, because as model inputs, these may give very different simulated output, compared with simulations based on the net change at a coarse scale (Kato et al., 2013;Wilkenskjeld et al., 2014;Fuchs et al., 2015).For example, a reported net increase in forest area of 10 km 2 may actually result from afforestation of 50 km 2 and deforestation of 40 km 2 .As input data to an ecosystem model, this might produce quite different results, compared to the parsimonious assumption (afforestation of 10 km 2 and no deforestation) (Levy & Milne, 2004;Krause et al., 2016).Over most of the globe, data on land-use change are typically limited in spatial and temporal resolution, and are typically represented by a time series of the area occupied by each land-use class (Rounsevell et al., 2006).Little information is available on the gross changes which bring about this time series (Prestele et al., 2017).The IPCC Good Practice Guidelines recommends the estimation of land-use change matrices for reporting GHG fluxes arising from land-use change (Penman et al., 2003).This provides explicit information on the areas which have changed from each land-use class to every other class.Whilst these matrices contain more information, they are only valid over the single time period for which they were derived, being a two-dimensional summary.For modelling over longer time periods, these are not very useful in themselves.
To properly represent the change in land use over time, we need a higher-dimensional data structure.
Land-use change is not easy to measure.A key problem is identifying change from repeated map or survey data, where the magnitude of the change signal is very small against the background noise of sampling and measurement error.Large censuses and careful survey techniques are required to distinguish true change from differences arising from measurement and sampling error (Fuller et al., 2003).A further problem is that information on land-use change at national scale typically comes from multiple disparate sources, which are often inconsistent with each other, using different land-use classifications and definitions (Phelps & Kaplan, 2017), arising from different thematic areas, and focus on different spatial and temporal domains, with different resolutions (Fisher et al., 2017).For example, land-use data in the UK are available from the agricultural census and surveys, the national forestry sector, the national mapping survey, as well as earth observation products such as Corine, MODIS and the CEH Land Cover Maps.However, no single data source provides a reliable estimate of land-use change with national coverage which extends suitably far back in time.A data assimilation approach is needed to make best use of the available data, so as to provide such a product.Existing methods ignore the large uncertainties which arise in estimating past land use change, and data assimilation approaches can explicitly address this issue.
In general terms, data assimilation is an approach for fusing observations with prior knowledge (e.g., mathematical representations of physical laws; model output) to obtain an estimate of the distribution of the true state of some phenomenon.It has become very commonly used in fields such as atmospheric and oceanographic modelling, and numerical weather prediction (e.g.Lunt et al., 2016).Various techniques are used, such as simulated annealing, ensemble Kalman filtering, and 4D variational assimilation.All of these can be seen as special cases within the Bayesian framework, where models, parameters and data are related in a formal way via Bayes Theorem (Wikle & Berliner, 2007).There are some significant differences in applying data assimilation in our land-use context, compared with atmospheric modelling.
Firstly, there is only a very simple model, compared with the complex physical models of the atmosphere or ocean.By contrast, the observational process by which the data are produced is extremely complex, compared with the simple observations of air or sea temperature or pressure.Also, we are predicting retrospectively (i.e."hind-casting") over many years in the past, rather than "nudging" forecasts as new data becomes available.
Our aim here was to develop a generic Bayesian approach, using multiple sources of data, to make spatially-and temporally-explicit estimates of land-use change.In a case study, we apply the approach to Scotland over the period 1969-2015.As an example application, we use a simple model of carbon fluxes following land-use change to show how uncertainties surrounding land-use change can be propagated through to model output.

Mathematical approach and notation
We represent land use u as a number of discrete states from the set {forest, crop, grassland, roughgrazing, urban, other} encoded as integers 1-6.106 At each time step, we have a square transition matrix . 112 At each time step, the net change in the area occupied by each land use is given by the gross 113 gains (the vector of column sums, G) minus the gross losses (the vector of row sums, L): where β ujt and i and j are the row and column indices.
We thus have three data structures, U, B, and A, which are inter-related by equations 1 -4.
U contains complete information about the system, which can be summarised in the form of A and B. B contains partial information about the system, which can be summarised in the form of A, but does not directly specify U.In itself, A does not directly specify either U or B, but can be used as a constraint in their estimation.
Multiple data sources are available which provide information in the form of these different data structures.Our approach here is to use equations 1 -4 as a simple model to relate the different observational data via Bayesian data assimilation in a two-stage process.Firstly, we

Data sources
We combined a number of data sources (Table 1) to describe the spatial and temporal change in land use in Scotland in the approach outlined above.A classification scheme was produced for each of these to aggregate the data into the broad classes used by classes (Penman et al., 2003).This was considered coarse enough that differences between 137 classifications could be aggregated into these six common classes, so that translation between

Data assimilation 144
Our data assimilation method proceeded as follows.1978, 1984, 1990, 2000, and 2007, and we interpolated linearly between survey years to produce an annual time series.We used the estimates derived in this way as our prior distribution of B. Each year, the mean of the prior distribution was taken to be the value of B from CS.The standard deviation σ of the prior distribution was estimated by applying a bootstrapping approach to the CS data (Scott, 2008).
• National Agricultural Census (AC) data provide annual records of the total area in the main agricultural land uses (Scottish Government, 2017).The Agricultural Census is conducted in June each year by the government agriculture department.Farmers declare the agricultural activity on their land in the form of ca. 150 items of data via a postal questionnaire.The results are collated at national scale.These are a long-running data set with near-complete coverage of agricultural land, relatively consistent over time, and are reported as national statistics and to the FAO.Hence it is desirable for our estimates of land-use change to be consistent with these data as far as possible.We therefore use these data as observations of A ut in the Bayesian framework, and predict ∆A ut from B t according to equation 4. The likelihood of the net change observed by Agricultural Census (∆A obs ut ) arising from normal distributions with means determined by equation 4 and the parameter matrix B is where ∆A pred A skewed normal distribution of this form (Azzalini, 2017) gives the likelihood of the gross changes observed as: where φ is the standard normal probability density function, Φ is the corresponding cumulative density function, and α is the skew parameter.Positive α produces a positive skew (when α = 0 we have the standard normal distribution).The parameter α can itself be estimated as part of the data assimilation procedure.
• Several data sources provide observations of U for one or more land uses at a restricted set of time points.We combine these into a single array U obs as follows.
-For an initial estimate of U, we use the Corine data sets for 1990, 2000, 2007, and 2012(European Environment Agency, 2016).For each grid cell, change between these years was assumed to occur at a random time within the interval, so that at national scale we effectively interpolate linearly.This produces U with complete UK coverage at annual resolution over the period 1990 to 2012.
-We overlay this with IACS data over the period 2004 to 2015 (Tomlinson et al., 2017) We can therefore add an additional term to the likelihood function which incorporates the comparison of the observations B obs with the values in the current parameter set B pred .
• To establish the posterior distribution, we use the Markov Chain Monte Carlo (MCMC) approach with the "DEz" algorithm implemented in the R package BayesianTools (Hartig et al., 2017).For each interval in the 46 year time series, an MCMC simulation was run, using the prior B t matrix from Countryside Survey, the observations of ∆A t , L t , G t for that year, and the observed B t matrix from Corine-IACS_NFEW.In practice, it is more convenient to use log-likelihoods, and our overall likelihood was the summation of log(L net ), log(L gross ) and log(L B ). Nine chains were used, with 100,000 interations in each.To establish the initial B parameter values for one of the chains, a least-squares fit with the ∆A was used.Other chains were over-dispersed by adding random variation to this best-fit parameter set.
• -EDINA Agricultural Census gives an estimate of ∆A at 2-km resolution.For each land use, an observed increase in area indicates the likely location of predicted gains.We therefore add a term to w which is proportional to ∆A.
-The CEH Land Cover Map (Rowland et al., 2017) gives an estimate of U t in 1990, 2000, 2007, and 2015 at high spatial resolution.Occurrence of a land use in the LCM suggests an area where gains would be more likely to occur.We add a term to w, based on occurrence of that land use in the LCM.
-Agricultural Land Capability Maps gives an estimate of how suitable land is for intensive agriculture, with a scale which ranges from good arable land, through intensive grassland and extensive grassland, to rough grazing.This scale can be translated into a probability of occurence for the land uses considered here, and added into the weighting of the sampling again.We use all the above information to produce many posterior realisations of U post , using the posterior B matrix and the sampling process described earlier.
Because the U data structure is large, we are limited in simulating many samples.It is therefore useful to summarise as the much smaller set of unique vectors and their corresponding areas.Our approach is to simulate 1000 samples, to calculate the unique vectors and their areas, and not to retain the larger data structure to reduce storage requirements.Another possible approach would be to simulate using only the MAP B matrix, and thereby generate the most likely realisations of U xyt , rather than the whole posterior distribution.

Carbon dynamics following land use change
We applied a simple empirical model of carbon fluxes following land use change, based on the UK LULUCF GHG inventory (Griffin et al., 2014).The soil component is based on the work of Bradley et al. (2005), and uses an analysis of the total soil carbon stock in a large number of soil cores, classified by land use and soil series.A linear mixed-effects model was applied to these data, to quantify the average effect of land use on soil carbon stock, treating soil series as a random effect.The model uses these mean values to represent the equilibrium soil carbon stock for each land-use class.When land use changes, the soil carbon stock moves towards the equilibrium soil carbon stock for the new land use.The soil carbon stock at location (x,y) and time t is given by: where C eq u is the equilibrium soil carbon stock for the current land use u, C xy,t−1 is the soil carbon stock at the previous time step, and k is a rate constant.The flux of carbon over the time step, ∆t, is given simply by difference: The above-ground component applies to the growth of biomass following afforestation, and uses the yield tables for British forestry produced by Edwards & Christie (1981), as interpolated and expanded to include non-merchantable timber biomass and wood products by Dewar & Cannell (1992).The mean change in above-ground biomass was assumed to be negligible in other land-use transitions in this simple model.

Results
Because of the availability of remotely-sensed data products, we are relatively confident in the present-day distribution of land use (Figure 2).This shows the concentration of urban areas in Scotland in the central belt, the restriction of cropland to the drier, flatter east coast, improved grassland mainly in the lowlands in the wetter south and west, and rough grazing and forestry sharing the Southern Uplands and Highlands in the north and west.
As an initial step in the data assimilationn process, a close least-squares fit to ∆A was achieved within a few tens of iterations, indicating that there were no particular numerical difficulties in estimating the B parameters.Standard measures were applied to assess whether the posterior distribution of B was suitably characterised by the output of the MCMC sampling.As well as inspection of the trace plots and the form of the distribution of the B parameters, we calculated the effective sample sample size, the acceptance rate, and various standard convergence diagnostics (Gelman & Rubin, 1992;Geweke, 1992;Raftery & Lewis, 1992).All of these showed satisfactory performance, that the MCMC chains converged, and that nine chains with 100,000 samples provides a reasonable estimate of the posterior distribution of B.     For cropland and improved grassland, CS and EAC show general agreement on the magnitude and pattern in area gained and lost to each land use (Figure 5 and Figure 6).An exception is an apparent anomaly in the early 2000s, when EAC gains and losses are both around 1000 km 2 higher than average for two years.This is not reflected in the net changes reported in the AC, so has to be treated with some caution.Reported gains and losses of rough grazing are much higher and very variable in EAC.This variability does not seem closely linked to the net change reported at national scale, so again, we treat this with some scepticism.There are no data on the gross gains and losses of urban and other land-use areas, as they are not covered by the AC or CS, and these terms are less well constrained.in the prior, the B matrices from CS which shows markedly higher crop to grass and grass to crop conversion rates over this time.
Figure 8 shows the 20 most frequent vectors more clearly, with each vector on a separate panel.This shows that 17 out of 20 involve transitions to or from rough grazing (which includes all semi-natural) land, which is the largest land use in Scotland by some way (around half the total area).Seven of these represent afforestation, which has mainly occurred on less productive, upland rough grazing land.Five vectors represent expansion of improved grassland on to rough grazing land.Vectors with two or more changes are less frequent, with none occurring in the top 20, but do represent a significant part of the total area (~8 % of the area undergoing change).In this simple soil model, land uses with higher equilibrium soil carbon than the average will tend to act as carbon sinks; those lower than the average will be sources.Carbon emissions from cropland increase as predominantly grassland is converted to cropland between 1970 and 1990.This then levels off as the cropland area remains stable or declines thereafter.Transitions to forest and rough grazing result in carbon sinks because they both have higher than average equilibrium soil carbon, and both show sizeable gross gains over the period.Rough grazing land also shows substantially larger gross area losses, but the associated carbon fluxes associated with this are attributed mainly

Discussion
The results show that we can provide improved estimates of past land-use change using multiple data sources in the Bayesian framework.The computation involved is quite feasible on a modern computer, requiring around three hours to estimate the parameters for a 46-year period.The output of the assimilation procedure provides vectors of land-use change in the form required for dynamic and process-based modelling, which we illustrate with the soil carbon modelling example.The main advantage of the approach is that it provides a coherent, generalised framework for combining multiple disparate sources of data.
As far as we are aware, there are no previous applications of formal data assimilation approaches to land-use change.However, some studies have addressed the same problem with related methods.Hurrt et al. (2006Hurrt et al. ( , 2011) ) used estimates of A together with estimates of wood harvest to predict B. The study was carried out at global scale at 0.5 degree resolution, and covered both historical and future scenarios for the period 1500-2100.To make the problem tractable, the transition matrix B was initially specified for only three land uses, so that a unique minimum solution could be found.Additional transitions associated with shifting cultivation and wood harvest were then calculated in a further step.They used a rule-based model which specified assumptions about the residence time of agricultural land, the priority of land for conversion to agriculture and for wood harvesting, and the spatial pattern of wood harvesting within a country.The distribution of land use over space and time U was not explicitly represented; instead, the area and age of "secondary" land in each grid cell was tracked in a book-keeping approach.However, because only a matrix is calculated at each time step, the approach does not produce explicit vectors of land use for dynamic modelling, and such things as rotational land use are not easily represented.Sensitivity to various assumptions was analysed, but the uncertainties associated with the input data and these model assumptions cannot readily be quantified.This approach yields essentially the same data structure as our method, and is wider in scope, covering all of Europe.
Our method represents an advance on this in several ways.Because the approach of Fuchs et al. ( 2013) is based on net change in areas at country scale, the extent of the true, gross changes will be under-estimated, possibly by orders of magnitude, and implicitly the B matrices are minimised.Our approach uses explicit observations of the annual transition matrices B as far as possible.Rather than regression relationships, our approach uses annual spatially explicit observations of where and when land-use change is likely to have occurred (based on CS, IACS and EAC).We use higher temporal and spatial resolution (annually, at 100 m) because this is possible with the data available in the UK, and with the limited spatial domain we attempt to cover.At continental and global scales, the same quantity and resolution of data is not available, and the computation issues become much larger.Our approach explicitly incorporates and propagates the uncertainty in the posterior distribution of B and predictions of A and subsequently modelled carbon fluxes.The uncertainty in land-use change is substantial, even in the UK where land management records are good.
Our methodology accounts for this uncertainty in a mathematically rigorous way (Van Oijen, 2017), and propagates this through to the subsequent modelling of other outputs, such as soil carbon fluxes.On a fundamental level, the Bayesian approach gives the correct theoretical answer to the data assimilation problem: if the observational error and prior are correctly specified and the posterior is adequately characterised by the MCMC sampling, then the posterior correctly represents the actual state of knowledge about the system parameters and predictions (Gelman et al., 2013;Reich, 2015).
We thus need to consider how well we can characterise the observational error, and the prior and posterior distributions.Establishing that the posterior distribution has been adequately characterised by the MCMC sampling is relatively straightforward.There are various criteria for assessing this (the effective sample size, and measures of MCMC chain convergence) which the results meet.In this study we chose to use an informative prior based on CS.This follows the way in which the data became available chronologically; these were the only data available with which we could estimate land-use change in the UK when an inventory of carbon emissions was first attempted (Cannell et al., 1999).The uncertainty in the prior distribution of B can be relatively well quantified, because considerable effort has gone into quantifying the likely level of error in the national-scale estimates of land use (Scott, 2008;Wood et al., 2017).The standard deviation σ of the prior distribution was most easily estimated by applying a bootstrapping approach to the CS data, but more advanced approaches have been investigated (Henrys et al., 2015).Alternative options for the prior are possible, and would be worth exploring further to examine sensitivity to the specification of the prior.Where little information is available, an uninformative prior is often used, either uniform, or exponentially declining to capture the parsimony principle that low values of B are more likely than high ones, all else being equal.More usefully, because we iterate over all years independently, we could form the prior distribution at time t from the posterior distribution for the previous year.In practice, we iterate backwards in time, so in fact the posterior at time t becomes the prior for time t − 1; this is mathematically simple but linguistically confusing.This approach means that information gained in the recent part of the time series is carried over into the earlier part of the time series.Subsequent estimates "borrow strength" from previous ones, in the Bayesian terminology.Currently, we do not use this approach because of the extra computation time this incurs, but methods to speed up this step can be explored.
Observational error can be difficult to estimate objectively and accurately, and often the σ terms are poorly known.Even in relative terms, it can be hard to judge the degree of certainty to place in different data sources, where observational error is not readily quantified.
In our case, we need to estimate the σ terms in the likelihood function (equations 5 -7) for the AC, EAC and IACS data.Spatial coverage in the data sets is similarly large so there is no clear a priori reason to trust one more than the other.However, there are reasons to prioritise the national-scale trends in AC over those from IACS, and to be cautious of the spatial patterns in EAC.AC is a long-established survey with relatively consistent methods, whereas IACS is a recent introduction, and the recording methodology has not been entirely stable over this period (for example, with changes to how much farm woodland is recorded).
It also attempts to collect a much higher level of detail (at the individual field scale), and this brings more potential for misclassification to appear as ostensible land-use change.However, with the limited information available, we cannot rule out that this is the more accurate data set, and that EAC and CS underestimate gross change.The accuracy of spatial information in EAC is limited by the way in which the data are collated, using postcodes of the land owner who completes the census return.Where large estates are owned, the correspondence between the centroid of the postcode district and the actual location of the land may not be very close.We therefore ascribe lowest uncertainty to AC, and higher but equal uncertainty to EAC and IACS data.In our Bayesian data assimilation procedure, IACS-based estimates of B are effectively down-weighted when they produce a mismatch with the national-scale AC trends.IACS coverage on forest, urban and other land is not large, and we would not expect accurate detection of changes in these land uses.
One of the main problems in land-use studies is that of classification.Depending on definitions used to delimit land-use classes, quite different areas may be calculated for the same nominal classes, and there is a real problem in combining data from different sources in that we may not be comparing like with like.Here, we minimise this problem by using a relatively coarse land-use classification, with only six classes.This would become more problematic if attempting to distinguish more refined classes.The computation time and difficulty increases with the square of the number of land-use classes, so there may be practical limits to the level of detail in the classification used, especially if applying on larger spatial domains.
An attractive feature of the Bayesian data assimilation approach is that additional data sources can be added to the process as they become available, without any major changes to software or step-changes in results.Several other data sources exist in the UK which could be incorporated.These include spatial data on the granting of woodland felling licenses, which would further constrain the likely location of deforestation, and national mapping agency data on urban expansion.As new satellite instruments come on-stream (e.g. from Sentinel and synthetic aperture radar), further remotely-sensed data products will become available which could be added into the estimation of A, B and U.In this study, we do not attempt to forecast future land-use change, but in principle this is simple with this methodology.If no new data are available, the posterior distribution will widen as future years are iterated over.
If scenario data were supplied, such as projected forest planting rates (G) or cropland areas required for food security (A), these could be used in the estimation of A, B and U in the same way as historical data.The method has applications in providing estimates of historical land use and land-use change input data for modelling work in many domains, including climate modelling (Lawrence et al., 2016), ecosystem and biogeochemical modelling (Ogle et

Figure 1 :5
Figure 1: Graphical depiction of a hypothetical 3-D cuboid U representing land use in space and time dimensions.Different colours show different land uses.
use a Bayesian approach to estimate the parameters in B, given prior information and partial observations of U and A. Secondly, we use the posterior distribution of B and spatial and probabilistic information on the location of land-use change to simulate posterior realisations of U. The maximum a posteriori probability (MAP, the mode of the posterior distribution) realisations represent our best estimate of land use and land-use change, given the available data.
ut is the prediction from equation 4 for the change in land use u at time t, and σ obs ut is the observational error in the Agricultural Census.So, we now have (i) a simple model which predicts net land-use change in terms of a parameter matrix; (ii) prior estimates of these parameters for each year from the Countryside Survey; and (iii) a function (equation 5) for the likelihood of the observations of net change given the model parameters.Combining these in Bayes Theorem, we can estimate the posterior distribution of the parameters, the transition matrix B. However before describing this, we can extend this simplest likelihood function by adding further sources of observational data.• The EDINA Agricultural Census (EAC) data (http://agcensus.edina.ac.uk/) provide additional information on land-use change, as they attempt to produce a spatially explicit version of the national-scale Agricultural Census data.Farm-level data is aggregated to 2-km grid cells, and data are available (or can be inferred) annually.While not containing explicit information on the actual land-use transitions, the resolution of the data is high enough that the net changes recorded each year in each 2-km cell may approximate the gross changes.In other words, because the data records the annual increases and decreases in land use across the grid of 2-km cells, the national totals of these increases and decreases gives an estimate of the gross change, the row and column sums of the transition matrix B, as well as the net change.When calculating the likelihood in our Bayesian framework, we can thus use the more informative observations of gross gains and losses (G and L) rather than just the observations of net change (∆A) from the national Agricultural Census.However, we know that the observations will tend to underestimate the gross change, because of the nature of the data reporting process: any counter-balancing gross change within the 2-km square is not included.To account for this, we can use a skewed normal distribution to represent this, such that predictions which overestimate the observations are more likely than underestimates.
Having established the posterior distribution of B, we use spatial and probabilistic information on the location of land-use change to simulate posterior realisations of U post .Starting with our best estimate of the near-present state of land use, U obs t=2015 , we work backwards in time.At each time step, we know the number of grid cells which need to change from land use i to land use j from the posterior matrix B t .For each i to j transition, we perform a weighted sampling operation to select this number of cells from those where U xyt = i.In choosing which cells to assign to j, we use the available data to calculate the probabilities which weight the sampling.Recall that U obs is given by the amalgamation of Corine, IACS and NFEW data.In the simplest case, the probabilities are determined only by this: all cells where U obs xyt = i and U obs xy,t−1 = j have equally high probability of being selected in the sample, and all cells where U obs xyt = i and U obs xy,t−1 = j have equally low (but non-zero) probability of being selected in the sample.This requires only a few simple rules to construct the probability weightings, w, for sampling cells for conversion from i to j: if U obs xy,t = i then w xy ← 0 else w xy ← 1 ∧ if U obs xy,t−1 = j then w xy ← 1 else w xy ← p mwhere p m is the probability of cells being misclassified in U obs , which we estimate to be 0.05.Sampling is done without replacement, so that a grid cell can only be selected once per year.To illustrate with an example, we start with our current map of land use, U obs t=2015 .Suppose our posterior estimate of B t determines that seven grid cells change from crop to grass, as we go back to 2014.Only cells which are crop in 2015 are valid candidates.Of these, those which were grass in 2014 (according to U obs ) will have high probability of being selected; others will have a low probability.If the posterior β post ijt area is lower than β obs ijt , not all the cells with high weightings from the above rules will be selected in the sample.If the posterior β post ijt area is higher than β obs ijt , additional cells, with low weightings from the above rules, will be selected in the sample.Thus, the cells which we are likely to change are those which are designated by U obs as crop in 2015 and grass in 2014.The effect of this is to generally recreate the spatial and temporal pattern seen in U obs (data from Corine, IACS and NFEW), but modified according to the extent of change estimated in the posterior B post .• As well as using the data from Corine, IACS and NFEW, we can also use other spatial data sets to inform the location of land-use change in our simulatations of the posterior U xyt .Any spatial data set which gives information on where and when a land use or land-use change occurs can be incorporated into the weighting used for sampling.Here, we used three additional data sets.

Figure 3
Figure 3 shows the Agricultural Census observations, and posterior predictions of the net change in area of each land-use class.The net change implied by the prior CS and IACS observations of B are also shown.The broad trends are: (i) an increase in forest cover due to sustained commercial forest planting; (ii) a corresponding decrease in rough grazing and semi-natural land due to expansion of forestry and improved grassland; (iii) an increase in cropland area between 1970 and 1990, with subsequent decline to the present day, due to changes in economic forces and subsidy incentives; (iv) an increase in grassland area since around 1990, partly corresponding to the reduction in crop area, and partly due to a general expansion on to rough grazing areas; and (v) a slow but consistent expansion of the urban area.These trends are picked up by the different sources of observations to some extent.The

Figure 2 :
Figure 2: Land use in Scotland in 2015 as estimated by the CEH Land Cover Map."Grass" comprises all improved and actively managed agricultural grassland."Rough" includes all rough grazing, unmanaged grassland and semi-natural land."Other" comprises barren areas such as montane and coastal areas.For legibility, we show this aggregated to 2-km squares, though the data are available at 250-m resolution 339

Figure 4 :Figure 5 :
Figure 4: Prior and posterior distributions of the transition matrix B, representing the gross area changing from the land use in each row i to the land use in each column j each year from 1969 to 2015.Red lines show the prior estimate from the Countryside Surveys.Pale blue points show estimates from IACS plus Corine and NFEW.The maximum a posteriori estimates after assimilating all data sources are shown in purple.The shaded band shows the 2.5 and 97.5 % quantiles of the posterior distribution.Note the y scale is different for each row.CS provided our prior estimate of B. Given the relatively small spatial coverage of CS, 340

Figure 6 :
Figure 6: Time series of the gross loss in area from each land use (A ut ) from 1969 to 2015, showing the observations, prior and posterior estimates.The shaded band shows the 2.5 and 97.5 % percentiles of the posterior distribution.

Figures 3 -
Figures 3 -6 show that there is considerable spread in the posterior distribution of B and predictions of ∆A.The 95 % credibility interval is typically of the order of 100 km 2 for the individual B parameters, and several hundred km 2 for the predictions of ∆A.The credibility intervals are smallest where multiple data sources agree on the nature of land-use change, and where the change is coherent across land uses.That is, an increase in one land use has to be balanced by a decrease in one or more other land uses.We have less confidence in predictions where the observed change in one land use is not compensated for by other land use changes.Credibility intervals in ∆A increase as we go back in time, because the uncertainty accumulates from year to year, although the increase has square root form rather than linear, Figure7and Figure8attempt to convey the detailed structure of the posterior U in a simple graphical summary.Figure7shows the 100 most frequent vectors of land-use change.Line thickness and opacity are proportional to the frequency (= area) of each vector, so that the dominant vectors are the most visually obvious.The plot shows that a wide range of land-use transitions occurs over the time period considered.Transitions from rough grazing to forest and to improved grassland are dominant.Bi-directional transitions between crop and improved grassland are particularly common in the 1980s.This comes from information

Figure 7 :
Figure 7: Trajectories of the 100 land-use vectors in the posterior U with the largest areas (excluding the six vectors which show no change).Each vector of land use is shown in a different colour, varied arbitrarily to differentiate different vectors.Line thickness and opacity are proportional to the frequency of (or total area occupied by) each vector, so that the dominant vectors are the most visually obvious.

Figure 8 :
Figure 8: Trajectories of the 20 land-use vectors in the posterior U with the largest areas (excluding the six vectors which show no change).Line thickness is proportional to the frequency of (or total area occupied by) the vector

Figure 9
Figure 9 shows the CO 2 flux resulting from land-use change over the 46-year period, derived from equations 8 -9 and the posterior distribution of U. The positive fluxes denote a gain to the terrestrial carbon stock, negative fluxes represent a loss to the atmosphere.We only represent land-use change from 1969 onwards here, but the effects on carbon flux are long-lasting.Hence, the carbon flux calculated here is initially small, and increases as the area having undergone land-use change accumulates over time.The accumulation of carbon in forest biomass (and wood products) following afforestation over this period is the largest term in these results.The forest planting rate has decreased markedly since 2005, giving the reduction in carbon sequestration in recent years.In this simple soil model, land uses with

Figure 9 :Figure 10 :
Figure 9: Net carbon flux from land-use change in Scotland over 1969-2015 showing the maximum a posteriori estimate and its 95 % credibility interval.The flux is attributed to change to each land-use class u.Positive fluxes denote a gain to the terrestrial carbon stock; negative fluxes represent a loss to the atmosphere.

Figure 11 :
Figure 11: Net carbon flux (in kg C m −2 ) from land use change in Scotland over 1969-2015 from the maximum a posteriori estimate of U .Positive fluxes denote a gain to the terrestrial carbon stock; negative fluxes represent a loss to the atmosphere.

Fuchs
et al. (2013) used a number of data sets, including that ofHurrt et al. (2006), to explicitly estimate the change in land use over space and time U for the whole of Europe at 1 km 2 resolution for each decade 1900-2010.Using logistic regression, they calculated "probability maps" for each land cover class, based on biogeophysical and socio-economic properties of each grid cell as explanatory variables for land use in 2000.For each decade and each country within the EU27, the net increase in the area of each land use (positive ∆A ut ) was allocated to the grid cells with the highest probability score for that land use.

can be derived from U t by comparison with the previous layer U t−1 . Each 110 element is given by 111
Bradley et al.

Table 1 :
Data sources assimilated in the estimation of land-use change in Scotland.
. The Integrated Administration and Control System (IACS) is a Europeanwide spatially explicit dataset at the field level that serves as a register of agricul-Again, this only has limited coverage, as it only covers forest land, but is given precedence in the case of conflict with the Corine/IACS data.We iterate over each time step to calculate B obs data are missing.Where there are conflicts with Corine, IACS data are given precedence because they are direct ground-based records.-Wethen add forestry data from the GB Forestry Commission (FC) National Forest Estate and Woodlands (https://www.forestry.gov.uk/datadownload), which records the location and planting date of forestry.t thus contains an observed estimate of the transition matrix for each year, from the combination of Corine, IACS and FC data.