A Linear Programming Model to Biophysically Assess Some Ecosystem Service Synergies and Trade-Offs in Two Irish Landscapes

Ecosystem service provisions are becoming more frequently used to assess land-use related conflicts in recent decades. This study investigates the current spatial and research information available to quantify ecosystem services relative to forest land-use planning in Ireland. A model is developed using the linear-programming method in Remsoft’s Woodstock platform. This model is applied to two case study areas in Ireland: Western Peatlands and Newmarket. Each case study area was chosen to assess a unique issue in the Irish and European context on the provision of ecosystem services. Western Peatlands was chosen to assess the effects of forest and alternative land-use options and Newmarket was chosen to investigate the effect of afforestation. The synergies and trade-offs of biophysically optimising the provisions of each ecosystem service are presented and discussed. The study quantitatively determines that trade-offs among provisions of some ecosystem services are required when optimising an ecosystem service while other ecosystem services are synergistic when the provision of a single ecosystem service is optimised.


Introduction
Deforestation, relating to land settlement and Ireland's colonial past, resulted in a forest cover of circa 1% around the beginning of the 20th century. The most recent National Forest Inventory was carried out in 2012 and estimated the total forest cover to have increased to approximately 10.5% [1]. Coillte (the Irish State Forestry Board, Newtownmountkennedy, Ireland), owns circa 53% of this, with the remainder owned by private landowners (ibid). Forests, public and private, have typically been planted on marginal land (Not suitable for agriculture but suitable for forestry) since the middle of the 20th century. As this marginal land is not suitable to grow high quality broadleaves, the majority of Irish forests comprise (exotic) conifers such as Sitka spruce (Picea sitchensis (Bong.) Carr), lodgepole pine (Pinus contorta) and Norway spruce (Picea abies (L.) H. Karst). These three species make up approximately two-thirds of the Irish forest estate [1]. However, many areas of Ireland have the biophysical potential for a range of land-uses and tree species.
Studies in a forestry context have suggested that proper sustainable forest management (SFM) will depend on the balance between its three pillars (environmental, social and economic) and their indicators [2][3][4]. Ecosystem services (ES) are a utilitarian concept which may prove useful to identify this balance [5]. A full history of the evolution of the ES concept is described by [6]. Briefly, the origin of the concept dates back to the late 1960s [7]. It was used in the 1970s to capture public interest in terms of biodiversity conservation [8]. Since then, the number of academic publications referring to the term ESs has gradually increased culminating in the Millennium Ecosystem Assessment; a major piece of work which provided an indication of the level of degradation of the world's ESs [9]. Quantifying ESs would make it possible to evaluate trade-offs and compatibility of ESs between different management scenarios in a given landscape. However, debate continues about the most appropriate method to define and quantify ESs with consistency. It is broadly accepted that localised studies will provide a more accurate assessment of ESs than on the global scale [10]. Managing for the provision of multiple ESs is one of the key challenges of ecosystem provision and state that trade-offs and synergies of ESs are particularly significant issues in addressing this at the landscape level [11]. DeGroot, et al. [12] have suggested indicators that link ESs to human well-being. In Ireland, political attention has focused on specific ESs, and land-use disputes have identified ESs that are strongly related to human welfare. Research in Ireland has investigated these [13][14][15][16] and most of the relevant ESs can be summarised into the categories outlined by DeGroot, Alkemade, Braat, Hein and Willemen [12].
Since the 1960s, many forests in Ireland have been established with a focus on generating revenue and this remains the main objective generally. This single financially oriented objective has guided forest decision support system (DSS) research in Ireland in the past. Nieuwenhuis and Williamson [17] developed a system for timber harvesting and sawmill delivery while minimising costs. More recently, SFM has come to the fore in Irish forest management and research effort has shifted. Forest research effort has become less focused on timber production with more of an emphasis now being placed on forests' non-timber benefits in parallel with timber supply, ESs. The potential to optimally group sub-compartments (the smallest unit of management used by Coillte) to create economically viable blocks for harvesting was explored by Nieuwenhuis and Tiernan [3]. They went on to implement some SFM constraints and compared their economic impact. PractiSFM, a simulation based DSS was developed to provide an approach for private forest owners to compile and forecast various management prescriptions while viewing the associated trade-offs of a range of ESs [2].
The objective of this study is to develop a DSS to quantify and evaluate the trade-offs of ES provisions at the landscape level over the course of a 70-year planning horizon, beginning in 2012. As forest rotations in Ireland are mostly 35 to 45 years, a 70-year planning horizon is sufficient to capture the full ES provision cycle. By using two Case Study Areas (CSAs) (i.e., the Western Peatlands and Newmarket) in which the developed model was applied, an indication of ES provision potential will be obtained given the specific biophysical conditions in each CSA, while excluding Irish forest policy, regulation and industry demand requirements. These provision levels will provide the production frontiers for each ES in each CSA, while the effect of optimising one ES on the provision levels of the other ESs will give an indication of their compatibility.
In order to investigate the long-term future of forest growth and the associated ES provision, planning on a strategic time scale is required. The linear programming (LP) component of Remsoft (Woodstock) (Remsoft, 2014.07, Fredericton, New Brunswick, Canada) was chosen as the most suitable tool to investigate the provision of ESs in the fragmented and complex landscapes chosen for this study. Optimisation will deliver the preferred management prescriptions available and will identify ES production frontiers, trade-offs and synergies associated with these prescriptions. Access to the Remsoft software in both the academic and industrial partners in the project facilitated the development of the model through the exchange of spatial and management information.

Western Peatlands
The Western Peatlands (WP) CSA is located in the northwest of Ireland (Figure 1). It is based on a Business Area Unit, a method of land division and management used by Coillte. There are eight Business Area Units in the Republic of Ireland. Their boundaries follow town land boundaries (as much as possible) and they were designed to be independently profitable units [18]. The location of the WP CSA means that wind exposure is a dominant factor for all land-use options. The east of the CSA is more fertile and has more potential for land-use applications than the west which enables the potential for a range of agricultural land-uses, i.e., dairying and beef production. However, afforestation on these better soils is a much smaller political issue in this CSA than the management of the forested peatlands. the WP CSA means that wind exposure is a dominant factor for all land-use options. The east of the CSA is more fertile and has more potential for land-use applications than the west which enables the potential for a range of agricultural land-uses, i.e., dairying and beef production. However, afforestation on these better soils is a much smaller political issue in this CSA than the management of the forested peatlands.

Newmarket
The Newmarket (NM) CSA takes in the north of county Cork (Figure 1). The central part of the CSA, known as the Duhallow region, is designated as a disadvantaged area vis-à-vis the allocation of EU structural funds. Dairying is the major farm enterprise in the NM CSA. Tillage is also quite substantial, but is confined to the most fertile sites in the eastern parts of the CSA. There are strong regional variations in terms of farm scales and incomes. Farms located at lower elevations that are mainly specialising in dairying and tillage production tend to be bigger and more economically viable than smaller farms located in the hills that rely on cattle and sheep farming. Forests are almost exclusively located in the less fertile, upland areas of the CSA. Considering that afforestation will be investigated in this CSA, both forestry and agricultural land is included. Statistics that provide a description of each CSA are provided in Table 1.

Newmarket
The Newmarket (NM) CSA takes in the north of county Cork (Figure 1). The central part of the CSA, known as the Duhallow region, is designated as a disadvantaged area vis-à-vis the allocation of EU structural funds. Dairying is the major farm enterprise in the NM CSA. Tillage is also quite substantial, but is confined to the most fertile sites in the eastern parts of the CSA. There are strong regional variations in terms of farm scales and incomes. Farms located at lower elevations that are mainly specialising in dairying and tillage production tend to be bigger and more economically viable than smaller farms located in the hills that rely on cattle and sheep farming. Forests are almost exclusively located in the less fertile, upland areas of the CSA. Considering that afforestation will be investigated in this CSA, both forestry and agricultural land is included. Statistics that provide a description of each CSA are provided in Table 1.

Decision Support System and Model
The model was developed in Remsoft, a Canadian developed software package. Woodstock [20] accommodates linear programming and uses an exogenous modelling method, meaning that spatial constraints are determined before the optimisation takes place. It is intended for long-term strategic planning, based on input in the form of ESRI shape files (ESRI, Redlands, CA, USA), ESRI geodatabases or database files. This methods section will focus on the data and data manipulation required to develop and implement the model in Woodstock's linear programming structure. Woodstock builds a matrix in a mathematical programming system, which is solved using the MOSEK mathematical programming solving software. This study used Windows (Microsoft, Redmond, Washington, WA, USA) 7 Professional 64 bit (SP1) operating system with an Intel ® Core™ i7-3930K CPU @ 3.20 GHz (6 cores with 2 threads per core) PC with 32 GB of RAM.

Coillte Collaboration
Coillte developed a Remsoft based, country-wide timber harvesting model. The objective of this model is to maximise Net Present Value (NPV) from the harvesting and delivery of timber products to mill gate, while encompassing access, ecological and environmental constraints in the form of financial disincentives. The financial objective of Coillte's model stems from its commercial mandate. The company provided this base model (in February 2012) for use in this study.

Remsoft Model Setup
There are three bodies of work required in order to produce a functioning model for Remsoft Woodstock ( Figure 2): Management prescriptions (Section 2.4); GIS setup and spatial setup in Remsoft (Section 2.5); and Ecosystem Service yields and forest inventory attributes (Section 2.6). Each element requires its own setup prior to implementation. The model accommodates a maximum of three forest species per polygon.

Remsoft Model Setup
There are three bodies of work required in order to produce a functioning model for Remsoft Woodstock ( Figure 2): Management prescriptions (Section 2.4); GIS setup and spatial setup in Remsoft (Section 2.5); and Ecosystem Service yields and forest inventory attributes (Section 2.6). Each element requires its own setup prior to implementation. The model accommodates a maximum of three forest species per polygon.

Management Prescriptions
Management prescriptions are choices available to the model. The range of choices available depends on site and stand conditions. These prescriptions include retention, site preparation, fertilising, forestation, thinning (including continuous cover forestry (CCF)) and clearfelling operations ( Table 2).

Management Prescriptions
Management prescriptions are choices available to the model. The range of choices available depends on site and stand conditions. These prescriptions include retention, site preparation, fertilising, forestation, thinning (including continuous cover forestry (CCF)) and clearfelling operations ( Table 2).

Retention
No anthropogenic intervention takes place. Forests are allowed to mature until the next year of the model.

Site Preparation
Mounding A requirement to prepare a site for establishment on infertile peats. Not prescribed to native woodland sites.

Fencing
An entire site was fenced for afforestation and half of this cost was assumed for fence maintenance at reforestation. This is based on area rather than perimeter.

Fertilising
Sites which have supported a previous yield class of 14 m 3¨h a´1¨year´1 or less for conifers, or 6 m 3¨h a´1¨an´1 or less for broadleaves, receive 2 fertilisation applications (one at establishment and one at year 5) to ensure establishment. Sites whose forest productivity needs to be determined (i.e., agricultural sites) were assessed for fertilisation based on a fertility model adapted from Farrelly, et al. [21].

Forestation
The grant premium categories for afforestation as outlined by the Irish Forest Service in 2012 were used for afforestation and reforestation (reforestation was not grant aided although the options provide a suitable variety of restocking options); with the addition of a reforestation option for WP. Species suitability was determined based on soil type and elevation. Any species with a moderate to optimal suitability for a particular soil type according to Horgan, et al. [22] was considered suitable to be planted under the soil type criterion. Species which had a tolerance of 3-5 according to Horgan, Keane, McCarthy, Lally and Thompson [22] were considered suitable for elevations up to 150 m for broadleaves and 200 m for conifers, while species tolerant of exposure (a rating of 2 or less) were also considered suitable for elevations up to 300 m. Only lodgepole pine (Pinus contorta) was considered suitable for elevations above 300 m for the purpose of reforestation. It is assumed that weevil and vegetation control is applied to all sites on which forest is established (afforestation or reforestation) and that these sites are manually inspected by a forester to ensure that establishment has been successful.

Native woodland
Species selection is based on a site's soil type as describe in Forest Service [23]; a site must be capable of growing trees without the requirement for fertilisation at establishment.

Harvesting
The species with the highest productivity dictates the thinning regime. Conifer stands must be suitably productive for thinning to be an option (i.e., YC ě 10). Conifers can begin thinning on a 5-year cycle (for a total of two or four thinnings) at the age of 16 or 22 years (the optimisation chooses the optimal). The thinning interval and potential start ages for broadleaves is YC specific. For a broadleaf species with a YC of 12, thinning takes place at 5 years intervals starting ages of 10 or 16, YC of 10: 10 years intervals starting at ages 12, 20 or 23, YC of 8: 15 years intervals starting at ages 20 or 25, YC of 6: 20 years intervals starting at ages 23, 28 or 33 and YC of 4: 25 years intervals starting at ages 30, 40 or 45. Thinning for CCF can begin at the age of 19 or 42 years (both ages are choices in for the optimisation) and thinnings take place at 5-year intervals for the duration of the planning period.

Clearfell
Clearfelling is an option when the species in a stand with the highest proportion of area is a conifer and exceeds 18 m in mean height or if the species is broadleaf and its age is ě 60 years.

Agricultural Land-Use
The agricultural land-use and the associated practices are maintained, i.e., ruminant or tillage production (Newmarket only).

GIS Setup
The spatial information is interpreted by Remsoft in a shapefile format. The spatial information and sources required to produce the shapefile for this study are presented in Table 3. The land-use datasets (i.e., the GIS vector layers that contain the Coillte forest, private forest and agricultural land) were merged in order to produce the land-use layer for this study. Then, GIS techniques and data related queries carried out in ArcGIS 10.1 (ESRI, Redlands, CA, USA) and Microsoft Access 2010 (Microsoft, Redmond, Washington, WA, USA), respectively, were used to assign the necessary attributes to each polygon in the land-use layer.

Quantifying Ecosystem Service Provision Levels
The origins of the values used in the model for growth, yield, habitat and sustainability parameters are described in this section, as well as their association with various management activities and states of the forest. Details of the metrics used to quantify ecosystem services are provided in Appendix (Tables A7-A13).
Potential end uses for timber produced in Ireland are sawlog, firewood, biomass for energy production, or the raw material for the production by Irish-based companies of a variety of panel products (i.e., medium density fibreboard, oriented strand board and moulded door panels). Experience in Ireland suggests that improper site selection, establishment and management will result in the large majority of Irish broadleaf stands not producing quality sawlog. For this reason, all broadleaf harvesting volume was assumed to be sold as firewood. The volumes are based on the static yield tables produced by the British Forestry Commission [24]. Spacing at establishment was fixed depending on species, and yields are dependent on species, yield class and thinning status. Each static yield table includes diameter at breast height (cm), stocking (stems¨ha´1) and harvesting volumes (m 3¨h a´1). Stand growth is quantified in m 3¨h a´1¨year´1.
There is more information available to quantify NPV (€) than any other ES and this is reflected in this study. This section outlines the revenue and costs of the forest and agricultural practices used in the model: Silvicultural establishment costs (Table A1): It is assumed that all stands have suitable access for timber extraction and transportation, therefore a road construction cost is not included. However, road maintenance, repair and upgrade costs are included and were assumed to be incurred at a rate of €1.50 per cubic metre harvested in clearfell or thinning.
Annual afforestation premiums (Table A2). Forest products the harvested volumes (Table A3): Divided into assortments using an assortment table based on average DBH ( [25] Revenue from the sale of standing timber is then calculated based on the constant market price for each product. Agricultural gross margin (Table A4): The marginal group was assumed to be located on farms designated as hill farms according to Teagasc's NFS definition while the gross margin for the profitable group was based on the gross margin for all farm productivity types, including hill farms. Hill farms were included as information was not available just based on the more productive farms. The net revenue from each year is then discounted to the beginning of 2012 at a rate of 5%.
There are several important species from an Irish political perspective, these are the native red squirrel (Sciurus vulgaris), hen harrier (Circus cyaneus) and deer (Dama dama, Cervus nippon and Cervus elaphus) require suitable foraging and cover habitat.
The conservation, enhancement or reduction in the provision of biodiversity ecosystem services was assessed based on species richness (SR) and habitat suitability. Scales from 0 to 10 were developed, where 10 is the most suitable habitat for a specific indicator species or biodiversity group. There are six biodiversity ES categories: the relative SR of nesting birds and SR of ground vegetation (Table A5), and the relative habitat suitability for red squirrel, hen harrier, deer foraging and deer cover. The habitat suitability scores were determined through a literature review.
Forest structure as defined by age was the primary stand level indicator of hen harrier habitat consistently referred to in the literature [26,27]. A distinction was made between first rotation and subsequent rotation forest in a circular published by the Forest Service [28] (Table A6). These ratings were verified by a hen harrier expert.
Red squirrel habitat suitability was determined based on the suitability of a tree species' seed as a food source, and the stand's development stage as an indicator of its ability to produce habitat and seed [29,30]. Bryce and Macdonald [31] present a list of suitable species for red squirrel habitat. This list was confirmed by Flaherty, et al. [32] and Waters and Lawton [33]. Species in the list were aggregated into the groups (Table A7) and the associated scores are presented in Table A8.
The research was mainly focused on forestry, with afforestation of agricultural land only considered in the Newmarket CSA. The matrix outlining the habitat suitability of non-forest land-uses was completed by a panel of experts, instead of being based on published results. The process involved three steps and was carried out for all six biodiversity ecosystem services (Table A9): 1.
Assign ratings out of 10 for permanence (i.e., frequency of likely disturbance) and botanical diversity that are unrelated to the ES scores.

2.
Group land-uses with similar characteristics into functional groups based on Step 1.

3.
Complete the biodiversity matrix. The scores were determined relatively to the habitat ratings of the forestry land-uses. If a non-forest land-use was considered to be more suitable for a given ES, all ratings were re-scaled so that the most suitable land-use type received a rating of 10. This scaling is reflected in all tables.
Research has indicated that agricultural land is only suitable for deer forage if there is adequate adjacent deer cover. Mysterud and Østbye [34] determined that only agricultural land within 200 m of a forest is suitable for deer forage. To account for this, a buffer was applied to all forest sites in the Newmarket landscape using GIS. The area of agriculture within this buffer was expressed as a proportion of the entire agricultural area (64% of agricultural land was within 200 m of a forest stand) and the deer forage habitat scores for agriculture were scaled proportionally. The ratings presented reflect this scaling (Table A10).
Phosphorus (P) concentrations [35] and nitrogen (N) concentrations have been recognised as key factors in freshwater eutrophication [36,37]. Suspended solid concentration is also an important factor in the assessment of water quality, as it can restrict sunlight from reaching photosynthetic plants, reduce dissolved oxygen [38] and, when sedimentation happens, can clog up spawning beds for salmonid fish stocks, which play a vital role in the freshwater pearl mussel's life cycle [39]. There was insufficient scientific information to quantify nitrogen or phosphorus release; however, it was possible to encompass a full range of information for water sedimentation risk. The revised universal soil loss equation was designed to estimate point source soil erosion [40]. Scheins (2001) adapted this equation for use with GIS techniques (Equation (1)). Finally, Sivertun and Prange [41] refined the techniques for estimating the equation's parameters and their model was used in this study. One aspect of the model they proposed was changed for this study. They categorised all agricultural, scrubland and bog with the same coefficient. Considering that there is a political requirement to establish watercourse buffer zones in Ireland to protect water quality, a separate coefficient was introduced for scrubland/bog/buffer zone, giving these land-uses a 20% lower sedimentation risk than undisturbed forest (Table A11). The water sedimentation risk in any one year (i.e., P in Equation (1)) is totalled for an entire CSA, area-weighted and divided by the maximum theoretical risk (i.e., assuming that the maximum possible risk in each polygon is present in either CSA). Water sedimentation risk is then expressed on a scale from 0 to 100. It is assumed that all forest management interventions are carried out in accordance with best forest practice.
where P = sediment loss risk, K = soil type, S = slope length factor, W = proximity to watercourse, and U = land-use.
Recreation scores for forest land were derived, on a scale from 0 to 10, from a study by Edwards, et al. [42], using the scores for the Great Britain region. The close-to-nature forests were assumed to be native woodland sites, the combined objective forestry (i.e., various management objectives can be combined in a manner that satisfies diverse needs better than through zoning [43]) was assumed to be CCF and the intensive even-aged forestry was assumed to be traditional rotation based plantation forestry (Table A12). Non-forest recreation scores were derived from these based on the relative permanence and structure of the land-use type (Table A13).
The carbon quantified in the model is that of whole standing trees. Whole tree biomass of softwood species was calculated based on average tree volumes, biomass expansion factors and root:shoot ratios developed by Levy, et al. [44]; the whole tree biomass is then multiplied by stems ha´1 to calculate stand level values. Aboveground biomass values for broadleaves were calculated on a stand level using a default equation developed by Teobaldelli, et al. [45] and belowground biomass (BGB) values were calculated using a default equation developed by Cairns, et al. [46], which estimates BGB based on aboveground stand level biomass estimates. Default basic density values were used from the intergovernmental panel for climate change [47]. A value of 50% was used to convert total biomass dry weight to carbon [48]. Carbon is quantified in millions of tonnes of carbon present in any given year (M T C). The difference between this stock from one year to the next is carbon flow and by optimising for carbon present over all years the cumulative flow is also optimised.

Area Scaling, Species Mixtures and Open Space
The ES scores presented assume that 100% of the area is covered with the land-use specified. This is typically the case for agriculture; however many forest stands require roads and rides for access and often have areas where forest growth has not been successful. To account for this, ES provisions were scaled in two steps. Initially, the stocking for some forest stands, particularly Coillte forest stands (but also private forests) has been determined through an on-site assessment to be lower than 100%. The second scaling method is based on the assumption that all forest stands are on average only stocked at 85%; this is common industry practice. The score for scrubland is assumed for the open space determined through both of these steps.

Biophysical Model and Scenario Objective
The biophysical model is applied to each ES with separate model runs to maximise and minimise the provision of each ES within the feasible solution space. A minimising scenario will find the smallest objective function value that is within this space. Whether the objective sense is maximising or minimising, the LP solver is always optimising the sum of the annual scores of the ES in the LP objective function over the entire planning period Equation (2). For reporting purposes, ES scores are then divided by the total area of the CSA to provide a relative area-weighted score. Each biophysical scenario was ran for 70 years (with 71 recording points, as ES provisions are reported at the beginning and end of the first year, i.e., January and December 2012). This length was chosen as a suitable compromise between computing power, strategic detail and sufficient time for forest growth and management actions to take place. There is just one constraint in this model. Its purpose is to ensure that the total area, specified as being present in all land-use categories at the beginning of the planning period, is accounted for in all future years of the optimisation (Equation (3)).
Max or Min ES x " Max or Min indicates whether the objective function sense is maximising or minimising; ES x = Ecosystem Service that the model can quantify; i and j are indexes for the year and development type (i.e., the smallest management unit used in the Remsoft Woodstock planning model) respectively; i = 1..y (the length of the planning period); j = 1, n (the number of development types); c ij = the objective function coefficient corresponding to the jth development type in period i, e.g., the hen harrier habitat score for a certain development type in period i; X j = the area (ha) of development type j; and D j = unaccounted area.

Results
The models produce a large quantity of output data, to account for this the results are presented in three forms which provide a good indication of the trade-offs associated with the optimisation of each ES individually: (1) the provisions of ESs at the beginning of the planning period; (2) the trends of ESs over the course of the planning period; and (3) a comparison of the average trade-offs between ESs.

The Provision of ESs at the Beginning of the Planning Period
In this section, a comparison will be made of the ES provision levels in the two CSAs in January 2012, before the commencement of any management interventions. This comparison only includes forested land before the optimisation begins, as agriculture was not included in the WP CSA. The ES provision levels are presented, on a hectare basis, to account for the WP CSA being larger in area than the NM CSA. Both CSAs have an age class structure typical of a rotation based silvicultural system where stands are clearfelled when they reach commercial maturity, and as a result ES provisions vary only subtly between CSAs. In this section, those ESs with significantly different provision levels between the CSAs are described, as well as those ESs that have unexpected levels of provision.
The WP's slightly more mature age class structure accounts for the greater provision of carbon and red squirrel habitat. More sunlight penetrates the mature forest canopy than the younger canopies of thicket stage or commercially maturing forests. The extra sunlight promotes ground vegetation growth on the forest's floor.
The area of open space is an unexpected factor that effects the provision of ESs. The higher hen harrier ES provision for the WP's more mature forest age class structure may be unexpected considering that the hen harrier ES scores are higher for young forests (Table 4). This higher level of ES provision is related to the area of open space and the difference in ownership structure, which together produce this unexpected result. The WP CSA has a 0.3% higher proportion of open space is a contributing factor towards hen harrier suitability as this type of land is much more suitable than forest land. In addition, it was assumed that all privately owned forests in both CSAs are on their first rotation, while state-owned forests were considered subsequent rotation forests. The hen harrier suitability rating of first rotation forest continues to be high for a longer period after establishment than for subsequent rotation forest according to the Forest Service [28]. The higher proportion of Coillte owned forest in the NM CSA means that more forest area is subsequent rotation forest. These factors produce a higher hen harrier habitat score for the WP CSA.
Recreation scores are higher in more mature forests, but even though the WP forests are more mature, the WP recreation score is lower than the NM score at the start of the planning period. Any forest that contains some proportions of conifer and broadleaf species is considered a mixed forest and has a higher score for recreation than a stand with just conifer species. Since afforestation subsidies became available for all privately owned land (before much of the afforestation of the NM CSA took place), a condition was that some diverse conifers or broadleaves must be planted with spruce or lodgepole pine monocultures. In contrast, when much of the WP was afforested (at an earlier stage under the Western Package scheme), and given that the WP is less suitable for broadleaves, conifer monocultures were established more frequently and, as a result, the WP CSA has a lower recreation score than NM.

Ecosystem Services Trends Throughout the Planning Period
The type of land-use that is typically chosen when optimising each ES is presented in this section. From this point on, the presented results include land that is forest and agriculture. Each CSA will be reported on separately rather than comparatively. If a maximisation ES is presented in this section, the minimising scenario will also be presented in the same figure and vice versa.
In the maxcarbon, maxgveg and maxrec scenarios, forest stands are retained past economic or silvicultural maturity. Although these retained stands enhance multiple ES provisions, they are not suited to providing hen harrier habitat (Figure 3). Optimisation of red squirrel habitat favours mature stands; however, it does not favour the forest species currently present in either CSA. Therefore, alternative forest species are chosen to develop optimal red squirrel landscape habitability. This is illustrated in Figure 4, where the hen harrier habitat, as a result of the establishment of Scots pine stands (through afforestation and reforestation) in the maxrsquirrel scenario, is provided at a moderate level for the first 25 years. When the Scots pine stands mature, they become unsuitable for the hen harrier. spruce or lodgepole pine monocultures. In contrast, when much of the WP was afforested (at an earlier stage under the Western Package scheme), and given that the WP is less suitable for broadleaves, conifer monocultures were established more frequently and, as a result, the WP CSA has a lower recreation score than NM.

Ecosystem Services Trends Throughout the Planning Period
The type of land-use that is typically chosen when optimising each ES is presented in this section. From this point on, the presented results include land that is forest and agriculture. Each CSA will be reported on separately rather than comparatively. If a maximisation ES is presented in this section, the minimising scenario will also be presented in the same figure and vice versa.
In the maxcarbon, maxgveg and maxrec scenarios, forest stands are retained past economic or silvicultural maturity. Although these retained stands enhance multiple ES provisions, they are not suited to providing hen harrier habitat (Figure 3). Optimisation of red squirrel habitat favours mature stands; however, it does not favour the forest species currently present in either CSA. Therefore, alternative forest species are chosen to develop optimal red squirrel landscape habitability. This is illustrated in Figure 4, where the hen harrier habitat, as a result of the establishment of Scots pine stands (through afforestation and reforestation) in the maxrsquirrel scenario, is provided at a moderate level for the first 25 years. When the Scots pine stands mature, they become unsuitable for the hen harrier.   In the maxtimber scenario, many stands are allowed to reach commercial maturity; in addition, all agricultural land and all clearfelled sites are forested with timber species that can produce volume most quickly, and the model chooses a commercial rotation-based system. Large clearfelling events are scheduled during the middle of the planning period in this scenario, increasing hen harrier habitat (Figures 3 and 4) but also water sedimentation risk ( Figure 5). Forests are then restocked and allowed to mature in order to optimise timber production through a large clearfelling event in the final year of the planning period. The maxbird, maxdeerc, maxnpv and maxdeerf ESs have a similar provision level of hen harrier habitat as the maxhh scenario, by providing habitat through the clearance of forest (Figure 4). The maxnpv scenario provides benefits to hen harrier habitat though clearfelling and foregoing reforestation. Deer foraging is provided by cleared forest, but more so by mature forest, while scores are particularly low for maturing forest. Therefore, in the maxdeerf scenario the model chooses younger forest stands to be clearfelled at the earliest allowable stage, while stands that are mature enough to benefit from the high ratings from the start of the planning period are retained.
In the NM CSA some agricultural land-uses have higher ES provision levels than any forest landuse. In the maxdeerf and maxbird scenarios, grassland is retained. The model chooses the higher habitat scores of forested land over some tillage land-uses for these scenarios (Figure 3). Intensive grassland- In the maxtimber scenario, many stands are allowed to reach commercial maturity; in addition, all agricultural land and all clearfelled sites are forested with timber species that can produce volume most quickly, and the model chooses a commercial rotation-based system. Large clearfelling events are scheduled during the middle of the planning period in this scenario, increasing hen harrier habitat (Figures 3 and 4) but also water sedimentation risk ( Figure 5). Forests are then restocked and allowed to mature in order to optimise timber production through a large clearfelling event in the final year of the planning period. In the maxtimber scenario, many stands are allowed to reach commercial maturity; in addition, all agricultural land and all clearfelled sites are forested with timber species that can produce volume most quickly, and the model chooses a commercial rotation-based system. Large clearfelling events are scheduled during the middle of the planning period in this scenario, increasing hen harrier habitat (Figures 3 and 4) but also water sedimentation risk ( Figure 5). Forests are then restocked and allowed to mature in order to optimise timber production through a large clearfelling event in the final year of the planning period. The maxbird, maxdeerc, maxnpv and maxdeerf ESs have a similar provision level of hen harrier habitat as the maxhh scenario, by providing habitat through the clearance of forest (Figure 4). The maxnpv scenario provides benefits to hen harrier habitat though clearfelling and foregoing reforestation. Deer foraging is provided by cleared forest, but more so by mature forest, while scores are particularly low for maturing forest. Therefore, in the maxdeerf scenario the model chooses younger forest stands to be clearfelled at the earliest allowable stage, while stands that are mature enough to benefit from the high ratings from the start of the planning period are retained.
In the NM CSA some agricultural land-uses have higher ES provision levels than any forest landuse. In the maxdeerf and maxbird scenarios, grassland is retained. The model chooses the higher habitat scores of forested land over some tillage land-uses for these scenarios (Figure 3). Intensive grassland- The maxbird, maxdeerc, maxnpv and maxdeerf ESs have a similar provision level of hen harrier habitat as the maxhh scenario, by providing habitat through the clearance of forest (Figure 4). The maxnpv scenario provides benefits to hen harrier habitat though clearfelling and foregoing reforestation. Deer foraging is provided by cleared forest, but more so by mature forest, while scores are particularly low for maturing forest. Therefore, in the maxdeerf scenario the model chooses younger forest stands to be clearfelled at the earliest allowable stage, while stands that are mature enough to benefit from the high ratings from the start of the planning period are retained.
In the NM CSA some agricultural land-uses have higher ES provision levels than any forest land-use. In the maxdeerf and maxbird scenarios, grassland is retained. The model chooses the higher habitat scores of forested land over some tillage land-uses for these scenarios (Figure 3). Intensive grassland-based agriculture is afforested at the beginning of the maxhh and maxdeerc scenarios in the NM CSA, which is reflected in an increase in hen harrier habitat. This hen harrier habitat provision gradually declines as these afforested stands (year 10 to 15 in Figure 3) mature. They are clearfelled as soon as they become eligible (period 30 to 35 in Figure 3). These stands are not reforested and provide very high hen harrier habitat scores associated with scrubland for the remainder of the planning period, i.e., afforestation is necessary as scrubland cannot be created straight from agricultural land.

Comparing the Trade-Offs between ESs
The intention in this section is to provide a representation of the ES provisions over the entire planning period. ES provisions can vary year to year and therefore the average level of each ES over the entire 70-year planning period is used. For each ES that is optimised, the average ES provision for that ES is improved compared to the 2012 score (Table 5). This indicates that the landscape composition in 2012 provides a mix of ESs rather than specifically targeting one individual ES. Furthermore, when comparing the average ES provision in each biophysical scenario with the corresponding 2012 level, it is clear that each CSA can facilitate increased provisions of all ESs included in this study. The optimisation scenarios and the resulting trade-offs between all ESs for both CSAs are presented in Table 6. This table depicts those ESs that have an average provision level within the more desirable half of the average biophysical range (i.e., within the more desirable half of the average maximum and minimum biophysical scenario results) when optimising each ES. For all scenarios, at least one ES provision level is outside of its desired range. This indicates that a reduction in some ES provision level is always the consequence of optimising other ESs. The number of ESs that remain within their desired range is reduced particularly when an ES requires unique habitat. The optimal type of habitat for red squirrel is mature Scots pine, which, in 2012 is a minor species in both CSAs. When Scots pine is planted for the maxrsquirrel scenario, the resulting provision levels for many of the other ESs provisions decline into the less desired range. This is also typically the case for timber production, carbon sequestration and recreation provision, i.e., a total of 6 or less ES provisions are in their desired range for both CSAs as a result of optimising these ESs. In the minh2o scenario, all agricultural land is afforested during the planning period and this makes more ESs less compatible than in the WP CSA, where all land is already forested at the beginning of the planning period. The young forests in NM are less suited to providing a broad mix of ESs. When the compatibility is assessed in the WP, recreation was found to be compatible with timber production. This is unexpected, as recreation requires mature forest while timber is produced when mature forests are harvested. At the beginning of the planning period, much of the WP CSA has mature forests, i.e., a higher recreation score than the afforested areas at the beginning of the planning period in NM. Recreation and timber production is therefore more compatible when scores are averaged over the course of the planning period in the WPs. Even though, for both CSAs, there is a large clearfell in the final year of the planning period resulting in a very low recreation ES score in the final year. Table 6. The trade-offs between ecosystem servicess when optimising each ecosystem service separately (minimising for water sedimenation risk and maximising for all others).

Discussion
The concepts of integrated forest, landscape and ES management are very new in Ireland and this study has introduced these concepts to a wide range of managers and stakeholders. The model described in this study is a first attempt to quantify a wide range of ESs simultaneously, based on Ireland's current climate and silvicultural systems. Up to now, research was done separately for different land-uses and no unified ES quantification system was available. The ES scales and scores developed in this study were produced within the academic domain, but it is expected that the model will be used in practical and operational settings in the coming years. Further development and refinement of the ES scales used in this study will benefit the approach and improve the current methods used to determine ES scales and scores.
The academic and industry stakeholders consulted at the local and national level have extensive knowledge of individual ESs quantified and of the CSAs investigated in this study. Based on the success of this study, the model will be used in a new international research project that involves many of the same stakeholders. In addition, Coillte, who also use Remsoft, is interested in incorporating several ESs into its management planning procedures.
It is important to note that the Irish model and outputs described in [49] are different from those presented in this manuscript. The Irish model in [49] involves a goal programming approach that incorporates constraints to reflect Irish forest policy and regulations and that aims to optimise stakeholder selected ecosystem provision levels. The model described in this manuscript is a more fundamental and pure biophysical model that is used to establish the absolute boundaries on the ecosystem service levels that the landscape can provide. These boundaries are then used as goals for the ecosystem services in the models used in [49].

Caveats
The following caveats to the method described should be noted: Spatial adjacency: The model developed in this study omits spatial adjacency constraints from within the optimisation algorithm. Without spatial adjacency, some polygons (or groups of polygons) may be too small/large, or too fragmented/homogenous, thus rendering the non-spatial solutions sub-optimal [50].
Fluctuations in ES provision levels: The scenarios produce the highest and lowest cumulative scores for each ES provision over 70 years. Within these cumulative scores, considerable fluctuations happen in the annual provision levels over the planning period. This may mean that optimal solutions may not necessarily sustain the populations living within these landscapes for the entire planning period.
Landscape heterogeneity: The methods used to quantify the ES scales and scores developed in this model were chosen because they were applicable to Ireland's landscape heterogeneity, climate and silvicultural systems. Much of the research used was based on data collected in the United Kingdom, (i.e., timber, carbon and some biodiversity) but also information from Sweden (i.e., water sedimentation risk). The ES scores used have been accepted by academic and industry experts, however the changes in heterogeneity that do happen may mean changes in the ES scores from those currently used.
Other forestry-related ESs: The ESs chosen for inclusion in this study were considered the most relevant to Irish forest policy development and forest management planning. There are of course other ESs that are influenced by forest management e.g., climate regulation, cultural heritage and identity and genepool protection. Their inclusion will be considered in future model development [12].
ESs not specific to forestry: In the NM CSA, agricultural land-uses were included in the model. However, several ESs whose provision levels are typically higher in agriculture than in forestry were not taken into account. Some examples of these include pollination, food production and the provision of fibre, fuel and other raw materials (other than wood) [12].

Conclusions
The results of this study show that in both CSAs there were ESs that were compatible with many others, while optimising some ESs resulted in large trade-offs with others. This result is consistent with studies in other geographical areas, which also investigated the simultaneous provision of ESs [3,51]. This study was carried out to determine the production frontier for each ES and the trade-offs with other ESs required to reach these provision levels. These biophysical ranges can assist in demonstrating to stakeholders the trade-offs and synergies that optimising one ES has on the other ESs quantified within a model. Such quantitative analysis capability will be useful when facilitating stakeholder-based management planning discussions. When combined with forest policy, land-owner behaviour and stakeholders' preferred ES provision levels, the findings of such stakeholder discussions will be beneficial towards developing landscape level integrated management and land-use plans. Author Contributions: Edwin Corrigan and Maarten Nieuwenhuis designed the study and interpreted the consolidated results. Edwin Corrigan wrote the paper. Maarten Nieuwenhuis contributed to writing the paper. Both authors contributed substantial data, meta-information, and result interpretation.

Conflicts of Interest:
The authors declare no conflict of interest.

Land-Use Category €¨ha´1¨year´1
Mixed A native woodland must be commercially mature, i.e., over 60 years old to receive the native woodland scores. Stand must be 'commercially mature'.   Table A11. Land-use coefficients for water sedimentation risk calculation.

Land-Use Co-Efficient Description
Agricultural land 0.010 Any land used for agriculture (grassland or tillage) Cleared forest 0.040 Forest sites cleared in the last 4 years, this includes restocked sites and the first two years of afforested sites Undisturbed forest 0.005 Any existing forest site that is 3 years or older Scrubland 0.004 Land not used for agriculture or forestry e.g., scrubland or bog