Impacts of land-use management on ecosystem services and biodiversity: an agent-based modelling approach

The science of ecosystem service (ES) mapping has become increasingly sophisticated over the past 20 years, and examples of successfully integrating ES into management decisions at national and sub-national scales have begun to emerge. However, increasing model sophistication and accuracy—and therefore complexity—may trade-off with ease of use and applicability to real-world decision-making contexts, so it is vital to incorporate the lessons learned from implementation efforts into new model development. Using successful implementation efforts for guidance, we developed an integrated ES modelling system to quantify several ecosystem services: forest timber production and carbon storage, water purification, pollination, and biodiversity. The system is designed to facilitate uptake of ES information into land-use decisions through three principal considerations: (1) using relatively straightforward models that can be readily deployed and interpreted without specialized expertise; (2) using an agent-based modelling framework to enable the incorporation of human decision-making directly within the model; and (3) integration among all ES models to simultaneously demonstrate the effects of a single land-use decision on multiple ES. We present an implementation of the model for a major watershed in Alberta, Canada, and highlight the system’s capabilities to assess a suite of ES under future management decisions, including forestry activities under two alternative timber harvest strategies, and through a scenario modelling analysis exploring different intensities of hypothetical agricultural expansion. By using a modular approach, the modelling system can be readily expanded to evaluate additional ecosystem services or management questions of interest in order to guide land-use decisions to achieve socioeconomic and environmental objectives.


Model and Landscape Set-up
This process defines all variables, brings in GIS data, and includes model-wide user controls. Table A includes descriptions of all general and landscape variables; model-specific variables are presented in the appropriate model section below. An image of the model interface is presented in Figure A. Predominantly native grasslands m 2 lc_120 Agriculture m 2 c-p Forest with >75% conifer cover, at least half of which is pine m 2 c-s Forest with >75% conifer cover, at least half of which is spruce m 2 d Forest with >75% broadleaf cover m 2 mx-p Forest where neither conifers nor broadleaf cover make up 75% of the forest. At least half of all conifer cover is pine. m 2 mx-s Forest where neither conifers nor broadleaf cover makes up 75% of the forest. At least half of all conifer cover is spruce m 2

ag-qs
Area of agricultural land in a cell m 2 for-qs Area of forest cutblocks in a cell m 2 hard-qs Area of hard linear features in a cell (e.g. paved roads) m 2

Forest Timber & Carbon Model
This model has been developed to simulate forest growth and carbon storage in Alberta's Green Area, as well as estimating the economic value of harvested timber and stored carbon. Critical input landcover data include forest stand type and age, areas ineligible for timber harvest (e.g. riparian buffers), the location of mills, and a cost surface representing the cost to transport timber from forest stands to mills for processing. The model operates at an annual time step, and each step involves several sequential processes.

Forest growth
Forest growth is based on standardized yield curves (as a function of stand age) developed for each of 5 forest types (spruce-dominated coniferous; pine-dominated coniferous; spruce-dominated mixedwood; pine-dominated mixedwood; and deciduous) in each of the 2 primary forested ecozones in Alberta (Boreal Plains and Montane Cordillera), for a total of 10 forest strata in the model. Growth curves were obtained from Natural Resources Canada (2006). The percentage of coniferous vs. deciduous tree volume in each stratum, as well as initial stand ages, were obtained from ABMI's enhanced vegetation layer, which in turn is based on Alberta Vegetation Inventory (AVI) and provincial forest fire records. Further, each forest cell is defined as either eligible or ineligible for harvest, based on standard criteria for Alberta (parks, other protected areas, and riparian buffers are ineligible).

Timber harvest
The study area is divided into Forest Management Agreement areas (FMAs) that are managed by different forestry companies, and the locations of mills present in Alberta today are included in the model. In doing so, the model accounts for the spatial and temporal patterns of timber production costs and revenue. Each mill is assigned an area of land) from which to harvest timber (i.e. the appropriate FMA), a production capacity level (based on this historical public record), a mill type (sawmill, pulp and paper, oriented strand board, or mixed), and a preferred timber type (conifer or deciduous). The harvest algorithm asks each mill to examine its allocated forest area and harvests stands, prioritizing the stands > 80 years of age that have the least cost distance to the mill (See Timber Cost Distance below). The annual allowable cut (AAC) for each FMA is allocated to its mills, and mills continue to harvest timber until the AAC is reached, less a user-defined percentage to account for the fact that actual wood harvest tends to be less than total AAC, based on historical records (AESRD, 2013). Once harvested, the stand age of each harvested cell is set to zero, and the adjusted harvested volume is sent to the mill.

Economic valuation of timber
Timber value is defined as the market value generated from a given volume of timber; that is, timber value is equivalent to the profits generated at the mill from selling timber products. Net profits for each mill are defined as mill revenue, minus timber delivery and operational wood processing costs. Mill revenue is determined based on the volume and type of wood products (including lumber, pulp and paper, oriented strand board, or veneer) typically produced at each mill. Provincial 1-year results for timber revenue are provided inTable B. The profits from each mill are then mapped back to the landscape, based on the amount of wood that flowed from each cell to the mill in a given time period.

Model Carbon Storage and Economic Value of Carbon
Forest carbon storage-versus-age curves were defined for each of the 10 forest strata listed above (5 stand types in each of 2 ecozones), based on a published forest carbon model (CBM-CFS3;Kull et al. 2011;Kurz et al. 2009). The process adopted was as follows: First, forest strata yield curves were entered in to the CBM-CFS3 model for each ecozone. Second, a stand level simulation is run for each stratum-ecozone combination. After a 200 year simulation, the model outputs a series of carbon stock data for each carbon pool (in tonnes per ha) for each stratum within each ecozone. From this output data, a series of "carbon yield" equations are developed by regressing carbon stocks against stand age (See example for two carbon pools in Figure B). These relationships are then coded in NetLogo so that carbon storage is tracked for each forest stand as they age. The model tracks 9 separate pools, which are subsequently grouped into the 5 pools used by the Intergovernmental Panel on Climate Change (IPCC) for reporting purposes ( Table C).
The economic value of stored carbon is based on a per-tonne price for CO 2 -equivalent, which is set by the user. The model calculates both the total value of stored carbon in forested areas, as well as the change in carbon storage over the model simulation period.

Figure B.
Example of carbon storage vs. stand age for two carbon pools, generated using the CMB-CFS3 model.

Timber Cost Distance
Each mill has its own cost distance raster, representing the estimated cost per unit volume to transport timber from each cell to the mill. All cost distance rasters are generated from a single cost surface, which represents the cost to transport timber across each cell (in $/m 3 ). The cost to transport timber is based on a fixed hourly rate of $130/hr, therefore the main factor is travel speed (Joerg Goetsch, pers comm.). The $130/hr rate coverts to a per-km cost for a truckload at different travel speeds; travel speeds were assumed for different types of roads (Table C).
For chips, there are 2.75 m 3 per BDT (bone dry tonne), 21 BDTs per truck, therefore 57 m 3 per tuck of hardwood chips (assume hardwoods are chipped at roadside). There is less volume per truckload of softwood, which is not chipped, so a correction factor of 1.623 is applied to the volume of softwood timber that is transported during this phase of the Forest Timber model.
For the cells that intersect road features, each cell was characterized by the length of the highest-speed road that occurs within it. For example, if a cell contains a paved road (speed = 90), it is assumed that the wood will travel along whatever the length of paved road lies within it, even if other lower-speed roads are present. Using the Cost Distance tool in ArcGIS 10.1, we generated a cost distance raster for each mill based on this cost surface (i.e. the cumulative cost of cells along the least-cost pathway from each forested cell to the mill in question). This value represents the monetary transportation cost to move timber from a cell to the mill for processing.

area-harvested-from-thispatch
Area of forest within a cell that is harvested ha

patch-haul-cost
Cost to transport timber from a cell to the mill that has harvested its timber $/m 3

is-harvested
True/false indicating if a cell is harvested True/false haul-cost-XXX Cost to transport timber from a cell to mill XXX (using "mill-namecode" to identify each mill) $/m 3

timber-esv
Calculates timber production value for each patch by multiplying npv-per-m 3 by my-patch-vol-harvested $ per cell

timber-evsp
Calculates potential timber production (i.e. unharvested eligible forest volume) value for each patch by multiplying npv-per-m 3 by total-patch-volume $ per cell

Carbon Model Variables aboveground-soft
Calculates above ground forest carbon from conifer (see appendix A)

Tonnes of Carbon aboveground-hard
Calculates above ground forest carbon from decid (see appendix A)

Tonnes of Carbon belowground-soft
Calculates below ground forest carbon from conifer (see appendix A)

Tonnes of Carbon belowground-hard
Calculates below ground forest carbon from decid (see appendix A)

Tonnes of Carbon dom-above-other
Calculates above ground carbon content of dead organic matter other than forest litter (see appendix A)

Tonnes of Carbon dom-above-litter
Calculates above ground carbon content of dead organic matter from forest liter (see appendix A)

Tonnes of Carbon dom-below-slow
Calculates below ground carbon content of slow decomposing dead organic matter (see appendix A)

Tonnes of Carbon dom-below-fast
Calculates below ground carbon content of fast

Pollination Model
The Pollination Model focuses on assessing the economic value added to canola production as a result of pollination by native bees. Using empirical relationships obtained from Morandin et al. (2006) in northwest Alberta, we parameterized a quarter-section-level model of canola yield as a function of bee abundance, which was in turn based on an empirical relationship with the amount of uncultivated land within neighbouring quarter-sections. Vector-based landcover and human footprint data obtained from ABMI were used to delineate the boundaries of crop fields and uncultivated areas. The advantage to using the vector-based ABMI dataset to delineate polygons is that they capture any small uncultivated parts of a field that can contribute pollinator nesting habitat that would not be recorded using a 30m resolution landcover raster.
In each time step of this annual model, the following processes are simulated:

Crop rotation
To obtain the location of canola fields, we used 4 years of annual crop maps (2009-2012; resolution 30m or 56m) developed by Agriculture and Agri-food Canada (AAFC 2012). These raster layers were used to determine the identity of crops within each agricultural polygon delineated by the ABMI vector landcover layer. In each year, the amount of each crop type in each cell changes based on these data, cycling through a 4-year rotation; after the 4 th year, the rotation starts again. where NL represents the area of natural and semi-natural land (measured in hectares) within the target cell and the 8 adjacent neighbouring cells. This includes all landcover types undisturbed by humans, including pasture but not annual cropland.

Calculate bee abundance and canola yield
The estimated bee abundance was then used to calculate the seed deficit SD as follows: = −12.54 + 1.29 0.48 Seed deficit is a proxy for yield, representing the difference between a canola pod's maximum possible seed set with the observed seed set. Based on data provided by Morandin & Winston (2006), each additional seed per pod represents an increase in canola yield of 73.3 kg/ha. The difference between the estimated seed set with pollinators present and the expected seed set if all pollinators were removed represents the fraction of canola yield attributable to insect pollinators.

Calculate pollinator value
The additional yield attributed to pollinators is multiplied by the market price of canola seed to obtain the monetary value of pollinators in a given model year. For multi-year simulations, the model calculates the net present value of pollinators based on a user-defined discount rate. Provincial modelled results for the portion of canola profits attributable to pollinators, based on four years of crop rotation data (2009)(2010)(2011)(2012) are provided in Table E.

Biodiversity Model
The biodiversity model estimates the effect of human land-use footprint on a variety of species. The biodiversity index, ranging from 0 -100%, represents the difference between reference (i.e. "defootprinted") and current abundance of each species, averaged across all species. The biodiversity model is a simplified version of ABMI's biodiversity intactness index (ABMI 2012) that is able to easily integrate with the other ecosystem service models.
ABMI's original biodiversity index analysis calculates intactness for each species as the ratio of current abundance to reference abundance (multiplied by 100%), or the inverse if the current abundance is greater than the reference abundance. Reference abundance is the predicted abundance of the species when the human footprint is set to 0 ("de-footprinted"). Current abundance is the predicted abundance at current human footprint levels. The same set of models is used for both predictions. The abundance metric for the biodiversity index is currently the probability of occurrence. Probability of occurrence is modelled as a binomial variable using a logit-link, with several trials per site (sample). Specifically, each site includes, 4 quadrats for plants, mosses, and lichens; 4 soil samples for mites; and 9 point counts for birds. Habitat elements are modelled as negative binomial counts with a log-link, a single-trial binomial with a logit-link for cover or log-normal, as appropriate for the different measures.
In addition to using footprint types as explanatory variables, the models also include covariates for vegetation types and geographic location. The habitat covariates for the boreal and foothills regions are currently vegetation types and ages of forest stands derived from ABMI's wall-to-wall vegetation mapping. Vegetation categories include pine, non-pine upland conifer (spruce), lowland conifer (black spruce/larch), deciduous, mixedwood, shrub, grass, non-treed wetlands, barren and water (the latter two are assumed to be 0). Age classes of forest stands are tracked by year and summarized into 20 year age classes, but combined into 3 broad groups in the modeling: 0-40yr, 41-100yr and 100+yr. In the grasslands and parklands, soil types are used as the habitat descriptors, because these are mapped even where agricultural footprint has replaced the native vegetation over large areas.
The current NetLogo model uses a simplified version of the species-specific intactness models, providing an estimate of how the overall biodiversity index (i.e. the average biodiversity index across all species) responds to footprint levels. We generated an equation to fit overall biodiversity index as a function of the percentage of different broad footprint types within a quarter-section. Models were developed for hard linear features (e.g. roads), soft linear features (e.g. seismic lines, pipelines, road margins), urban/industrial areas, agriculture, and forestry cutblocks. Each equation is of the form Quarter-section biodiversity index = 100 + a*pcHF + b*ln(pcHF+1) + c*pcHF 2 where pcHF is the percent (0-100) of the human footprint type in a quarter-section, and a, b and c are model coefficients specific to each footprint type. These fit the observed results very accurately, but they apply to individual footprint types. In quarter-sections containing multiple footprint types, the equation is used sequentially, from the most-severe footprint type to the least-severe type (in descending order of severity, hard linear, soft linear, urban/industrial, agriculture, and forestry), in the following procedure: Based on a random sample of 10% (n = 11007) of all cells in the North Saskatchewan Watershed region, the simplified model fits the full model very well (r 2 = 0.94).

Water Purification Model
The water purification model represents selected hydrological processes in Alberta watersheds related to non-point source pollution through surface flow and erosion. In particular, the model was designed to identify source areas of pollutants, important areas for removing pollutants, and impacts to downstream water users. Note that this model only includes surface flow and run-off, and does not incorporate ground water flow. This dynamic model involves several processes:

Model Set-up
This step brings in all the relevant spatial data, including landcover, precipitation, river network, and monitoring points of interest. All figures provided in this document depict the North Saskatchewan and Battle River Watersheds in Alberta, Canada.

Precipitation and Surface Flow
A rainfall event is simulated, where the average annual amount of precipitation falls on each cell and runs off downslope to the river network, eventually flowing downstream to the river outlet.

Pollutant Loading and Deposition
Pollutant loading is based on nutrient export coefficients (N, P, TSS) or RUSLE (sediment) combined with annual precipitation. Deposition is based on landcover types.

Stream Flow & Tracking Loads
The flow of surface water, nutrients, and sediments through the river network is modelled to calculate the cumulative annual load of each pollutant throughout the river network.

Source and Sink Areas
After overland and stream flow are complete, nutrients and sediment that reach the river network can be traced back to their origin points on the landscape.

Model Set-up
A river network for a major watershed was created using the flow accumulation tools in the Spatial Analyst (Hydrology) toolset in ArcGIS, based on a DEM of Alberta ( Figure D). The resulting network was visually compared to a GIS layer of streams in Alberta to assess its accuracy; any incorrect river connections caused by importing a linear river layer at an 800m resolution were manually corrected. This layer is imported into NetLogo such that any 800m cell that intersects a river is classified as a "river" cell. A network is created so that each river cell tracks its upstream and downstream neighbours, allowing for water and information to travel throughout the network. Finally, water monitoring stations are created to facilitate tracking of water flow and quality metrics at points of interest, such as mouths of tributaries, inlets of important lakes, river outlets, and municipal water treatment plants.

Figure D.
North Saskatchewan and Battle River Watersheds and river networks in Alberta, Canada.

Precipitation and Surface Flow
Precipitation is based on mean annual precipitation across Alberta, and is represented as a set of "raindrop" water agents with a specified volume of water, corresponding to mean annual precipitation falling on each patch of land. Runoff for each patch is calculated as a percentage of total precipitation, based on runoff coefficients for each landcover type (Table H; Donahue et al. 2013). Because each cell contains multiple landcover types, each cell's runoff coefficient is calculated as the area-weighted average of runoff coefficients for each landcover category. After the precipitation event, the volume of water that runs off each patch of land moves to the adjacent patch with the lowest elevation, and this downslope movement continues until each raindrop agent reaches a river cell.

Pollutant loading and deposition
Nutrients (nitrogen, phosphorus, and total suspended solids -TSS) and sediment are loaded into surface water flow through two different modelling processes, outlined below.

Water Quality
Nutrients (nitrogen and phosphorus) and TSS are loaded into water flow based on export coefficients (measured in kg•ha -1 •mm -1 of annual precipitation) calculated for major landcover and human footprint types (Donahue 2013). Each 800m cell in in the model includes multiple landcover and human footprint types, so pollutant loading is calculated as an area-weighted average across each landcover type present within a cell. Thus, the amount of pollutant k initially loaded into raindrop agent j (P jk , measured in kg) at origin cell c is calculated as:

= ∑ • •
Where L is the number of unique landcover types in the origin cell of the raindrop j, E kl is the export coefficient for pollutant k in landcover type l (measured in kg•ha -1 •mm -1 of annual precipitation), A l is the area of landcover type l within a cell (measured in ha), and W is annual precipitation (measured in mm).
Nutrient removal occurs during overland flow, where a percentage of each raindrop agent's pollutant load is removed as it flows across downslope cells before reaching the river network. Pollutant removal percentages are assigned for each landcover type, and the total amount removed by a cell is the areaweighted average of the landcover types comprising the cell. Thus, the amount of pollutant k in raindrop agent j after it has passed through adjacent cell c + 1 is calculated as: Where D l is the percentage of each pollutant removed from landcover type l. This process is repeated as the raindrop agent passes through all downslope cells until it reaches the river network. Example maps are provided for phosphorus loading ( Figure E) and removal ( Figure F).

Sediment
Sediment erosion into surface water flow is calculated using the Revised Universal Soil Loss Equation (RUSLE;Renard et al. 1997). Sediment erosion from each patch is calculated (in tonnes) as

= • • •
Where R is the rainfall erosivity factor, K is the soil erodibiity factor, LS is the slope factor, and C is the area-weighted cover factor that accounts for the influence of landcover types (e.g. vegetation, human use, etc.) on erosion. The R value was set at 0.350, following guidance from Wall et al. (2002). The K factor was based on soil texture categories (Wall et al. 2002); soil information was obtained from the Soil Landscapes of Canada version 3.2 (REF). The LS factor was calculated based on a digital elevation model using the Terrain Analysis tools in SAGA GIS (Bohner & Selige 2006). The C factor was based on generalized C values for Alberta generated by Wall et al. (2002). We calculated the area-weighted Cvalue for each cell as

= ∑ •
Where C l is the cover factor for landcover type l, A l is the area of landcover type l within a cell, and L is the number of distinct landcover types within a cell. In addition to these terms, RUSLE typically includes a term for the support practice factor ( P), which represents practices undertaken to modify runoff such as changing grade, flow pattern, or the direction of runoff. Common examples of support practices include cross-slope cultivation, contour farming, stripcropping, and terracing (Wall et al. 2002). P values represent a percentage reduction in erosion, and range from ~0.1 (90% reduction in erosion) to 1.0 (no change) if no support practices are employed. Given the lack of spatially explicit data on support practices in Alberta, this term was omitted in the present model (corresponding to a P factor of 1.0).
An example output map of sediment generation is depicted ( Figure G). Sediment deposition occurs via the same process as removal of water pollutants.

Stream Flow & Tracking Loads
Once raindrop agents have reached a river cell, all of the water flow and quality information they contain (flow, nutrient load, and sediment load) are transferred to the river network. Because each node in the river network is connected to its upstream and downstream neighbours, water flow and loading information is transferred downstream, such that each cell in the river network is able to track the total, cumulative amount of flow and loading that passes through that point ( Figure H). Monitoring stations that have been set up along the river network track the same metrics; this allows for easy export of data from these points of interest.

Source and Sink Areas
Once surface and stream flow are complete, several metrics are calculated for each cell in the study region. Each raindrop agent keeps track of its origin cell and its route through other downslope patches until it reaches the river network. Therefore, the final amount of pollutant load in each raindrop, after it has completed its flow to the river, is mapped back to its origin cell, providing a map of actual pollutant supply from each cell ( Figure I). Combining the phosphorus supply and retention maps provides an estimate of the net contribution of a given cell to the phosphorus reaching the river network. Combining the supply map with the retention map ( Figure F

Water Model Calibration
We calibrated the water model based on cumulative annual phosphorus loading data obtained for 8 points in the headwaters of the North Saskatchewan watershed. We conducted a global calibration by varying the export coefficients; we obtained 10 parameter sets using Latin hypercube sampling, and the parameter set that performed best for each natural region. After this procedure, the Spearman's Rank correlation between the calibrated model results and observed phosphorus loads at the 8 monitoring points was 0.33. Amount of preciptiation that infiltrates into soil mm channel-elevation Channel elevation data from channel-dataset applied to patches. A more refined elevation that better captures channel features m rainfall Annual precipitation mm k-factor soil erodibility used in RUSLE n/a ls-factor slope length factor used in RUSLE n/a c-factor cover factor used in RUSLE n/a hydro-elevation Modified elevation data derived by averaging elevation and channel-elevation (to represent elevation of water column) m slope slope calculated from elevation data % runoff-coefficient percent of rainfall that leaves a cell via surface runoff % hillslope-sedimentgenerated eroded sediment generated at each cell calculated using RUSLE Volume of water currently on a cell m 3 load-n Total nitrogen released in surface runoff on a cell kg load-p Total phosphorus released in surface runoff on a cell kg load-tss Total tss released in surface runoff on a cell kg load-n-perha Nitrogen released as runoff per hectare kg/ha load-p-perha Phosphorus released as runoff per hectare kg/ha load-tss-perha TSS released as runoff per hectare kg/ha n-deposited Nitrogen originating from upslope that is retained by a cell during overland flow kg p-deposited Phosphorus originating from upslope that is retained by a cell during overland flow kg tss-deposited TSS originating from upslope that is retained by a cell during overland flow kg sediments-deposited Sediments originating from upslope that is retained by a cell during overland flow tonnes n-supply Nitrogen originating from a cell that eventually reaches the kg