Protecting wildlife habitat in managed forest landscapes—How can network connectivity models help?

Industrial forestry in boreal regions increases fragmentation and may decrease the viability of some wildlife populations, particularly the woodland caribou, Rangifer tarandus caribou. Caribou protection often calls for changes in forestry practices, which may increase the cost and reduce the available timber supply. We present a linear programming model that assesses the trade‐off between habitat protection and harvesting objectives by combining harvest scheduling and optimal habitat connectivity problems. We formulate the habitat connectivity model as a network flow problem that maximizes the amount of habitat connected over a desired time span in a forested landscape, while the forestry objective maximizes net undiscounted revenues from timber harvest subject to even harvest flow and environmental sustainability constraints. We applied the approach to explore the trade‐off between caribou habitat protection and harvesting goals in the Armstrong‐Whitesand Forest, Ontario, Canada, a boreal forest area with prime caribou habitat. Our model also incorporates Dynamic Caribou Harvesting Scheduling (DCHS), a harvest policy currently in a place in Ontario that aims to balance the forest management and caribou protection goals in northern boreal regions. In our study area, the implementation of DCHS appears to have relatively minor impact on timber supply cost. By comparison, maximizing the protection of caribou habitat would lead to a noticeable increase of the mill gate timber cost by $3.3 m−3 on average, while enabling habitat protection in an additional 5.0%–9.5% of the range area. Our model is generalizable and can be adapted for assessing habitat recovery and harvest goals in other regions.

protection goals in northern boreal regions. In our study area, the implementation of DCHS appears to have relatively minor impact on timber supply cost. By comparison, maximizing the protection of caribou habitat would lead to a noticeable increase of the mill gate timber cost by $3.3 m −3 on average, while enabling habitat protection in an additional 5.0%-9.5% of the range area. Our model is generalizable and can be adapted for assessing habitat recovery and harvest goals in other regions.

Recommendations for Resource Managers:
• Incorporating the concept of long-term habitat connectivity into forest planning can help reduce the negative impacts of harvest activities on caribou populations.
• Prioritizing habitat connectivity leads to a small increase in the overall harvest area because harvest has to be allocated to less productive and more geographically isolated sites to protect prime wildlife habitat containing old conifer stands.
• Maximizing the habitat protection would lead to a noticeable increase of the timber supply cost (by $3.3 m −3 on average), while enabling moderate increase of the protected habitat area (i.e., an additional 5.0%-9.5% of the range area).
• Implementation of Dynamic Caribou Harvest Schedules, which is the current harvesting policy in Ontario's boreal forests when caribou populations are present, causes only a minor increase of the timber supply cost in our study area. undisturbed forest areas. In particular, woodland caribou populations (Rangifer tarandus caribou) have been declining in areas of industrial forestry operations (Brandt, Flannigan, Maynard, Thompson, & Volney, 2013;Venier et al., 2014). Harvesting creates a network of clear-cuts in early successional stages, which increases the abundance of competing deer and moose populations and attracts predators of caribou (James & Stuart-Smith, 2000;Latham, Latham, McCutchen, & Boutin, 2011;Wittmer, Sinclair, & McLellan, 2005). Woodland caribou is a threatened species in Canada (SARA, 2002) and declining numbers in its populations pose a serious conservation problem (Festa-Bianchet, Ray, Boutin, Cote, & Gunn, 2011;Hebblewhite, 2017;Hebblewhite & Fortin, 2017). Recovery efforts for protecting caribou populations aim to create larger regions with undisturbed habitat and eliminate open spaces that serve as movement corridors for predators (EC, 2011(EC, , 2012ECCC, 2017).
Protection of woodland caribou habitat is a long-term policy that focuses on minimizing human activities that cause fragmentation, such as timber harvesting. However, protecting wildlife habitat may reduce the area of forest available for harvest and increase the cost of timber supply. To plan effectively, decision-makers must be able to assess the interplay between the extent of harvest activities and measures to protect caribou habitat (Martin, Richards, & Gunn, 2016;McKenney, Nippers, Racey, & Davis, 1997;McKenney, Mussell, & Fox, 2004;Felton et al., 2017;Ruppert et al., 2016). Ideally, the habitat protection measures should have minimal impacts on forestry activities, but as a practical matter, caribou populations and harvesting co-occur across many parts of boreal Canada.
Recently, we proposed a network-based approach to solve a habitat connectivity problem for woodland caribou (Yemshanov et al., 2019). Under this approach, a fragmented forest landscape is depicted as a network of habitat patches (nodes) interconnected by arcs indicating potential movement corridors for animals. Connectivity between adjacent habitats is formulated as a network flow problem through this partially connected habitat network. The work follows from Sessions (1992), who proposed the formulation of the connected habitat problem as a Steiner network. It also relates to Jafari and Hearne (2013) and Jafari, Nuse, Moore, Dilkina, and Hepinstall-Cymerman (2017), who proposed a network flow problem to establish a contiguous protected reserve, as well as patches m to patch n if patch n has suitable habitat present over a continuous time span, T min or longer. We depict the connectivity between patches m and n as a bidirectional pair of arcs, mn and nm, which indicate potential movement of caribou individuals between m and n in both directions. Caribou movement through the network of patches N is conceptualized as a flow through a subnetwork of connected nodes over a continuous timespan T min or longer. We introduce the nonnegative variables y nm and y mn to characterize the bidirectional flow through arcs nm and mn connecting a pair of adjacent nodes n and m with habitat.
Caribou require suitable habitat to support their foraging and reproductive behaviour. Each forest patch n may have the amount of suitable habitat. The level of habitat suitability to support caribou individuals depends on the composition and age of the forest in n (Ferguson & Elkie, 2004a, 2004bJohnson, Parker, Heard, & Gillingham, 2002;OMNRF, 2015). Clear-cut harvesting temporarily destroys the suitable caribou habitat as it reduces the amount of local foraging resources. Early successional stages following the harvest attract deer (Odocoileus spp.) and moose (Alces alces L.) populations followed by predators (black bears (Ursus americanus Pallas) and wolves (Canis lupus L.), which further increases the predator pressure on caribou (James, Boutin, Hebert, & Rippin, 2004;Latham et al., 2011;Wittmer et al., 2005). A harvested patch is expected to become suitable again for caribou in 40-60 years, as forest stands mature and adequate vegetation cover is restored. For the current work, we assume that decisions to harvest forest stands in patch n are binary with no partial harvesting options. Harvest decision implies clearcutting all stands in n and resetting their forest age to 0 after the harvest. We assume a patch will undergo a natural regeneration after the harvest.
For each patch n, we define a set of harvest prescriptions i, i = 1, …, I, where each prescription defines a possible sequence of harvest events over a planning time horizon T, including a scenario without harvest. We enumerate all possible prescriptions that can be assigned to forest patch n by a set of binary vectors of length T, p ni = {(1,0, …, 0), (0, 1,…0), …}, p ∈ P. The elements of each vector denote the harvest or no harvest conditions in a particular time period t, t = 1, …, T. A binary variable x ni , x ni ∈ {0,1} selects whether a patch n follows a C a n a d a U S A Forested areas Caribou ranges, no industrial forestry Caribou ranges, industrial forestry Industrial forestry, no caribou harvest prescription i with a vector of harvest times p ni . Only one harvest prescription can be selected for a patch.
The time since harvest and tree species composition defines the amount of suitable habitat in patch n in period t as well as the continuous time span t, t + 1, t + 2, … during which the habitat remains suitable. For each prescription i, the parameter b nit defines the amount of suitable habitat that could support caribou individuals in patch n in period t, t ∈ T (see the description of the parameter b nit in the Data section). A binary parameter, λ nit , identifies the presence or absence of suitable habitat in n in prescription i in period t (λ nit = 1 when b nit > 0 and forest stands in n are older than 40 years, and λ nit = 0 otherwise). We use the λ nit values to estimate the parameter τ ni , τ ni ∈ [0; T], which defines the maximum number of consecutive time periods when the habitat remains suitable over a planning horizon T in patch n in prescription i. The habitat capacity in patch n in prescription i in period t is defined as follows:

| Habitat connectivity problem
Our habitat model tracks connectivity between patches with suitable habitat that remains connected for continuous time span T min or longer, T min < T. A binary variable w nm , w nm ∈ {0,1}, defines whether a species flow is established between patches n and m throughout the continuous timespan T min or longer. We find a network of connected patches in the landscape that maximizes the total habitat amount in the network. To ensure their spatial connectivity, we formulate a network flow problem that injects flow into the network of selected patches. Each selected patch must receive flow from at least one neighbouring selected patch. We introduce a root node, n = 0, that is used to inject the flow into the habitat network. The root source node 0 is connected to all other nodes (patches) 1, …, N and can inject the flow into any node in the network (Figure 2). Adding a root node is a standard approach in formulating network flow models that aim to find a spatially contiguous set of land units (see Jafari & Hearne, 2013).
We assume that connected nodes can be fed by no more than a single incoming arc from one node, and the flow only comes to a node from one source. This assumption prevents the creation of loops in the connected network. A node n with incoming positive flow from nodes m (i.e., with one of the w mn values set to 1) becomes a part of the connected network, that is, In Equation (2), set N n defines nodes m = 1, …, N which are connected to node n and can transmit flow to n, including the root node {0}. Since only nodes with suitable habitat available over a continuous timespan T min can be connected, node selection must include the selection of a harvest prescription i (which defines the forest stand age and so the node's habitat suitability status). We define a binary variable z ni , z ni ∈ {0,1}, as the product of the selection of flow to a node n and harvest , so z ni = 1 only if a prescription i is selected for node n and there is an incoming flow to n. Constraints (3)-(5) linearize the product of two binary variables as: We formulate the habitat connectivity problem as selecting a subset(s) of fully connected nodes to maximize the total amount of habitat that is connected over T min continuous periods (or longer) in a landscape N over a planning horizon T, that is, where the penalty decision variable P 1 , P 1 ≥ 0, defines the number of nodes connected to the root node 0 above one node, and f 1 is the scaling factor, subject to constraints (3)-(5) and: Node 0 (source of the flow) Arcs nm connecting adjacent nodes n,m Nodes n (forest patches) Arcs 0nused for injecting flow from root node 0 to a sub-graph of connected nodes Selected connected nodes Injecting flow to a sub-graph of connected nodes F I G U R E 2 The network flow model concept. Arrows show the universe of arcs, or potential connections between the neighboring patches in the habitat network. Light-shaded arrows show connections from the root node 0 to nodes n in the network, which are used to inject the flow into the network. Bold arrows in red show the flow injected into the habitat network; thin arrows in red show the flow through the selected connected nodes. Patches in red show the selected connected nodes For a selected node n, constraint (7) ensures that the amount of flow coming to the node is equal to the amount of flow outgoing from the node plus the fulfilled demand capacity at n. The demand capacity is 1 if a node is selected and zero otherwise. Term N n + defines the set of nodes m that are connected to node n and can receive flow from n. Constraint (7) ensures the flow balance through a node n and its connectivity with the other nodes. Constraint (8) helps avoid cyclic connections and ensures that the flow only comes to a node n from at most one source. Constraint (9) defines the penalty variable, which specifies the number of connections to the root node 0 exceeding one. Maximizing the objective value (6) with the penalty P 1 in place minimizes the number of connections receiving the flow from node 0 above one. The penalty in the objective function is required because each connection from the root node 0 creates a separate contiguous network. Compared to a formulation that sets a fixed number of connections to the root node (e.g., one), the penalty formulation makes it easier to find feasible solutions in fragmented landscapes where creation of a single contiguous habitat network is impossible. The coefficient f 1 defines the magnitude of the penalty term in the objective function equation. Setting the f 1 value sufficiently high yields a solution with a single fully connected network, after first finding a feasible solution with multiple connections to the root node.
Constraints (10) and (11) ensure agreement between the flow selection variable y nm that defines the amount of flow through an arc nm and the arc selection binary variable w nm . Constraint (10) ensures that the flow between nodes n and m is zero when the arc nm is not selected. Constraint (11) ensures that the arc selection variable w nm is zero unless a positive flow is established from node n to node m. Constraint (12) ensures that the flow to a node n can only be established if the node retains suitable habitat over a continuous timespan T min or longer. The parameter τ ni in Equation (12) defines the longest time span when suitable habitat is available in site n in prescription i and was derived as follows. Consider the following prescription example in site n: 10-year time periods, t: 1 2 3 4 5 6 7 8 9 10 Harvest events --1 -------Habitat availability, λ nit : 1 1 0 0 0 0 1 1 1 1.
In this prescription, the longest continuous time span the habitat remained suitable over T periods is four periods, that is τ ni = 4. The τ ni values were calculated when generating the universe of harvest prescriptions I. Constraint (12) also relates the binary flow selection variable to patch n, w mn to the timespan parameter τ ni and allows the connections to node n if the node is assigned the harvest prescription i with the τ ni value equal or above the minimum time span T min .

| Harvest scheduling problem
Forest patches (nodes) in the network can be harvested for timber. The allocation of harvest maximizes the net revenue, subject to a target volume of harvested timber in each period t, even harvest flow constraint in consecutive periods t and t + 1 and a constraint that maintains a minimum average age of forest stands in the area N at the end of the planning horizon. We adopt the harvest scheduling Model I formulation (Johnson & Scheurman, 1977;M. E. McDill & Braze, 2000;M. McDill et al., 2002M. McDill et al., , 2016Martin, Ruppert, Gunn, & Martell, 2017). The model considers an area of N forest patches over a planning horizon of T periods. For each patch n, a set of possible harvest prescriptions i, i ∈ I, defines the sequences of all harvest actions over T periods including a no-harvest scenario. As defined above, binary variable x ni selects the harvest prescription i for a patch n. We only consider clear-cut harvest, which is the most common harvest type in boreal Canada (NFD, 2019). A forest stand can be harvested after it reaches a minimum harvest age of k years or older. Each patch includes only one stand characterized by age, tree species composition and a forested area, a n , that could be harvested for timber. For patch n, prescription i defines a sequence of harvest times with volumes of harvested timber V nit and a precomputed net revenue R ni from harvesting that timber over period T minus harvest, hauling and mandatory postharvest regeneration costs.
We define Q t as the volume of timber harvested in area N in period t and Q t min and Q t max as lower and upper bounds on the harvest volume in period t. We also define ρ n as the unit price of timber harvested from a patch n net of harvest and hauling costs. Even harvest flow over consecutive planning periods t and t + 1 is enforced by a proportion ε that limits the change of the harvest volume proportion in consecutive periods t and t + 1 by 1 ± ε. We also add a minimum bound for the average age of forest stands in the managed area at the end of the planning horizon T, E Tmin , and define E ni as the forest stand age in a patch n at the end of the planning horizon if prescription i is applied. The harvesting problem is defined as maximizing the net revenues, R ni , associated with managing the forest over T periods, that is, The harvest revenue R ni is calculated as the undiscounted net cash flow associated with harvesting forest in patch n over the entire planning horizon T in prescription i. In our formulation, the harvested volume is constrained by even harvest flow and target harvest volume constraints, so maximizing revenues minimizes the per-unit cost of harvesting a target volume of timber, subject to constraints (14)-(17). The revenue R ni is calculated as the timber value net of harvest, hauling and postharvest regeneration costs, e n : For each site n, the R ni values were precomputed for every prescription i before optimization. Since our goal was to assess the impacts of long-term habitat protection, the use of undiscounted cash flows enabled similar handling of harvest allocation and habitat protection over the long term. The use of discounting would prioritize short-term harvest revenues over long-term cash flow, and so would misrepresent the long-term environmental sustainability of the harvest patterns and their costs.
Constraint (14) ensures that each patch with harvestable forest is assigned just one prescription. Note that the set of prescriptions I includes a no-harvest scenario with zero revenues. Constraint (15) ensures that the harvest volume for each period stays within a range [Q t min ;Q t max ]. Constraint (16) specifies that the harvest volumes in consecutive periods t and t + 1 do not deviate by more than upper and lower bounds 1 ± ε. Constraint (17) sets the average age of forest stands at the end of the planning horizon T to be greater or equal to the minimum age target E Tmin . This constraint ensures that a portion of the old-growth forest stands is left unharvested, which prevents overharvesting.

| Incorporating the Dynamic Caribou Habitat Scheduling (DCHS)
In 1996, the Ontario Ministry of Natural Resources and Forestry (OMNRF) adopted a new sustainable harvest approach in areas where caribou populations are present (FPAC, 2018;OMNR, 2009;Racey et al., 1999). The DCHS is a part of forest management plans approved by the Province of Ontario. For a 10-year planning period, a DCHS concentrates harvest areas and subsequently aggregates clear-cut disturbances, while minimizing densities of logging roads. Aggregating the harvest in a few relatively large regions helps to avoid creating a pattern of small cut blocks across large areas.
During the planning process, the selection of DCHS regions is done for the entire harvest planning horizon. Periodically, a forest plan is updated after major disturbances when the current plan becomes no longer useful for harvesting or wildlife protection (Baskent & Keles, 2005;Martell, Gunn, & Weintraub, 1998). Currently, DCHS regions and their harvest times are delineated statically during the planning process and do not incorporate dynamic scheduling, which may be needed in response to changing landscape composition or timber demands. Regardless, aggregating the harvest within a DCHS region helps maintain connectivity between undisturbed habitats outside the harvested region but also offers some flexibility for reallocating cut blocks within the region.
We incorporated the DCHS principle into the harvest scheduling model by way of two scenarios. Our first scenario uses a static arrangement of DCHS regions and restricts the harvest schedule to fixed time periods and DCHS regions according to the provincial DCHS plan (Figure 3). Based on data received from OMNRF, we defined a set of static DCHS regions s, s ∈ S. Each region has a specified timeline (a 10-year period) when stands in that region can be harvested. We assume that fine-scale allocation of harvest within a particular DCHS region would maximize the net harvest revenues, as described above in the harvest planning problem (13). For each patch n and time period t, a binary parameter, δ nt , defines whether a patch n can be harvested in period t according to a static DCHS plan (δ nt = 1) and δ nt = 0 otherwise. Constraint (19) forces the harvest in period t to regions as prescribed by the provincial DCHS plan, that is,

| Dynamic allocation of harvesting in DCHS regions
Static delineation of DCHS regions and harvest schedules may not be sufficient to guarantee long-term connectivity of caribou habitat. One reason for implementing DCHS is to ensure that the regions harvested in a given time period are surrounded by an adequate area of undisturbed habitat. In short, DCHS aims to avoid situations when harvest occurs in adjacent DCHS regions at the same time. We capture this aspect with our dynamic DCHS scenario, which only uses the geographic boundaries of the DCHS regions but allows dynamic selection of harvestable DCHS regions over time. We introduce a set of constraints to penalize the allocation of harvest in adjacent DCHS regions in the same period t.
We introduce a binary decision variable, l st , to define harvest in DCHS region s in period t. The l st values are estimated via constraints (20) and (21), that is, where θ ns is a binary parameter that indicates whether a site n belongs to a DCHS region s (θ ns = 1 and θ ns = 0 otherwise) and M is a large integer value. For each time period t, we then define an S × S adjacency matrix of binary parameters, α sq , where α sq = 1 if the neighboring DCHS regions s and q share a common boundary and α sq = 0 otherwise. We then introduce a binary decision variable L sqt , s, q ∈ S, t ∈ T, to define the occurrence of harvest in a pair of neighboring regions s and q in period t (i.e., L sqt = 1), and L sqt = 0 otherwise. The L sqt value is a product of decision variables l st and l qt , which indicate harvest in regions s and q in period t, and is linearized via constraints (22)-(24): We then define a penalty decision variable associated with scheduling harvest in adjacent DCHS regions in the same period, P 2 , P 2 ≥ 0 via constraint (25), that is, Constraints (20)-(25) create a harvest pattern that minimizes harvesting in adjacent DCHS regions in the same period.
A related goal of DCHS is to concentrate the harvest in a particular period in a few large regions while keeping the other regions intact. Thus, we need to ensure that the harvest is concentrated in as few DCHS regions as possible while avoiding harvest in adjacent regions in the same period. We add a decision variable P 3 , P 3 ≥ 0, to the objective function equation to penalize the total number of DCHS regions harvested in the same time period t above one, that is, We also need to define a minimum threshold area, ϕ min , that can be harvested in DCHS region s if a decision is made to harvest that region in period t. We define the range of possible harvest areas in DCHS region s as {0, [ϕ min ; ϕ s ]}, where ϕ s is the total harvestable area in DCHS region s. This condition is established by introducing a binary decision variable, d st , and a disjunction constraint (27), that is, where μ nit is a binary parameter, μ nit ∈ {0,1}, that indicates the presence of harvest in patch n in prescription i in period t (μ nit = 1) and μ nit = 0 otherwise. The μ nit values are generated from the harvested volume data (i.e., μ nit = 1 when V nit > 0 and μ nit = 0 otherwise). The model also requires a masking constraint to ensure that no harvest occurs outside of the designated area, that is, where ξ is a binary mask indicating that a patch n belongs to a harvestable area (ξ = 1) and ξ = 0 otherwise.

| Combining the harvest and habitat connectivity objectives
We combine the harvesting and habitat connectivity problems in a single objective via relative weights. Our full objective maximizes the weighted sum of the connected habitat amount in area N and the net harvest revenues over T planning periods, that is, A normalizing factor γ rescales the harvest revenue term in Equation (29) to make it roughly the same order of magnitude as the amount of connected habitat. Scaling coefficients f 1 -f 3 adjust the relative weights of penalties P 1 -P 3 in the objective function equation. The scaling factor f 1 should be high enough to push the model to create a single (or minimum possible number of) contiguous subgraph(s) with the connected habitat.
Solving the objective functions (29)-(31) with different objective weights F enables exploration of the trade-off between maximizing the amount of protected habitat versus maximizing harvest revenues. The solutions for end-points of this trade-off when F values are set to 0 or close to 1 depict the most distinct policies when habitat protection or harvest revenues are prioritized, which we have further explored in our study. We composed the model in the General Algebraic Modeling System (GAMS, 2019) and solved it with the GUROBI linear programming solver (GUROBI, 2019). Table 1 lists the model parameters and variables. We ran the model on a HP Gen 10 workstation with dual Xeon Gold processors for 72 hr or until reaching a 0.5% optimality gap (whichever came first).

| Case study
We applied the model to examine harvesting and caribou protection strategies in the Whitesand-Armstrong Forest Management Unit (FMU) in northwestern Ontario, Canada (Figure 1). The area is adjacent to Wabikimi Provincial Park and includes portions of both the Nipigon and Brightsand caribou ranges (CPAWS, 2009;OMNR, 2012). The area has been moderately fragmented by logging, with timber delivered to mills in Thunder Bay, Ontario. Currently, a plan to build wood pellet and cogeneration plants in Armstrong, Ontario is under consideration by the provincial and federal governments as well as the Whitesand First Nation (Bieler, Trush, & Jakob, 2019). If established, a wood pellet plant and a cogen facility would ensure energy independence of local communities and create a sustainable market for fiber. However, large-scale harvest operations may increase forest fragmentation and cause the decline of local caribou populations. Protection of sensitive caribou habitat has been proposed as a management tool to help prevent further decline of caribou populations in the region (Neegan Burnside Ltd., 2014) but has to compete with forestry activities. Figure 4 depicts key spatial inputs including habitat intactness (Figure 4a), timber volume in current conditions (Figure 4b), hauling cost (Figure 4c,d), harvestable area and suitable habitat amounts (Figure 4f-h) and the distribution of high-use (category 1) and seasonal-use caribou habitat (category 2, Figure 4i). Lower and upper bounds on harvest volume over a period t Q t min , Q t max ≥ 0 a n Forest area in a node n a n ≥ 0 V nit Volume of merchantable timber available for the harvest at a node n in period t in harvest prescription i

| Data
We divided the study area landscape into 1 × 1 km patches (Table 2). For each patch, we estimated the potential amounts of suitable caribou habitat b nit for each harvest prescription and forest age by combining two approaches. First, we modelled the current habitat distribution in the study area ( Figure 4g)  . Suitable habitat amounts were calculated as a function of current forest composition, recent disturbances, esker and anthropogenic linear feature density and land cover type using resource-type selection coefficients developed by Hornseth and Rempel (2015) from the location data. Models were updated and extended to the lower boreal region, including the Brightsand Caribou Range (Rempel & Hornseth, 2018). The future distribution of suitable caribou habitat was based on a boreal caribou habitat model for Ontario's Northwest Region (Elkie et al., 2018), from which we predicted amounts of suitable habitat based on a combination of land cover composition and forest age (Figure 4h and Table 3). For every 10-year forest age class, for each habitat type (such as useable, preferred, and refuge habitats), if present in patch n in period t, a score of 1 was assigned, and a total habitat suitability value was estimated as the sum of these scores. When forest patches included a mix of different land cover types, the total habitat suitability value was estimated as a weighted average of scores for individual cover types and their corresponding areas. We assumed that forest stands regain suitable habitat status 40 years after harvest. We then estimated the habitat suitability values as a weighted average of the suitability values based on the current habitat distribution and the future land cover and age composition. The weighting factor for estimates based on the current distribution was set to 1 in time period 1, and linearly decreased over the planning horizon to 0 in period T. The sum of the weighting factors for the estimates based on the current distribution and future land cover/age was set to 1.
The study area may also experience other anthropogenic disturbances that are undesirable for caribou populations. We adjusted the suitable habitat amounts b nit using a habitat intactness coefficient that accounted for human-mediated disturbances in the area of interest. For each patch n, we calculated habitat intactness values for all combinations of harvest prescriptions i and time periods t. We estimated intactness by averaging the three criteria which negatively affect the amount of suitable habitat: the area proportion of nonlinear anthropogenic  Table 2). Previous assessments of long-term caribou movement patterns in northern Ontario suggested that caribou tend to select habitat at broad scales (i.e., in the 5,000-10,000 ha range) rather than finer scales (Hornseth & Rempel, 2015). Therefore, we estimated the presence of T A B L E 2 Summary of data assumptions

Assumption Description
Data spatial resolution (sites n) 1 × 1 k m suitable habitat at coarse scales using a grid of 5,000-ha hexagon cells and then interpolated the habitat values to a grid of 1 × 1-km map cells. For the current caribou distribution, the habitat selection and resource utilization were modelled using the approach presented in Hornseth and Rempel (2015). This methodology is consistent with the general habitat description for caribou (OMNRF, 2015). Similarly, the suitable habitat estimates based on future land cover composition and age were aggregated from the spatial resolution of individual land cover polygons in Forest Resource Inventory data at the same 5,000-ha resolution (Figure 4g). We used road network data from the CanVec database (NRCan, 2019) to estimate hauling costs, assuming an on-site harvest cost of $15 m −3 and delivery of hardwood timber to Armstrong, Ontario (the closest market, which has a proposed pellet plant) and softwood to the Resolute Forest Products mill in Thunder Bay, Ontario (Figure 4d,e). The hauling costs were based on typical estimates for northern Ontario conditions (Maure, 2013) and included the delivery cost with a hauling rate of $90 hr −1 , assuming a 40 m 3 truckload, 1 hr waiting time and an overhead cost of $4 m −3 ( Table 2). The cost of accessing remote sites was incorporated by increasing the per-unit hauling cost value from sites without road access in direct proportion to the distance to the nearest road. Based on discussions with specialists from the Ontario Ministry of Natural Resources and Forests, we adjusted the hauling cost value so the area-wide timber cost per unit is approximately 25% higher than the unit cost value without accounting for road access.
The starting values for stand age, merchantable timber volume and land cover composition were estimated from Ontario's Forest Resource Inventory database (OMNRF, 2018). To estimate the future timber volume in the harvest prescriptions, we used yield curves for northwestern Ontario from a recent timber supply study (McKenney et al., 2016). These yield values were adjusted by the projected annual losses of forested area due to fire disturbances using fire regime zones from Boulanger, Gauthier, and Burton (2014). The minimum age of harvest k was set to 70 years. We assumed that the area-wide mean forest age at the end of the planning horizon must be equal to or greater than the current mean age (i.e., approx. 80 years). We set the even harvest flow bounds to ±2% and the harvest planning horizon T to 100 years with 10 × 10-year time steps. Table 2 summarizes key data assumptions.

| Forest management and habitat protection scenarios
We evaluated six harvest and caribou protection policies for a set of harvest targets between 0.05 and 0.3 M m 3 -yr −1 . Two management policies were represented by groups of "harvest priority" and "habitat protection priority" scenarios. "Harvest priority" scenarios maximize the net harvest revenues without prioritizing the protection of caribou habitat and have the scaling factor F in the objective function set to zero. "Habitat protection priority" scenarios set the scaling factor to 0.99 and maximize the amount of caribou habitat that can be kept connected over T min periods or longer, while giving low priority (0.01) to maximizing the revenues from harvest. Both scenario groups meet the harvest target [Q t min ;Q t max ]. Each scenario group included "no DCHS," "static DCHS," and "dynamic DCHS" scenarios for problems 1-3, as defined earlier (Figure 3). We also examined optimal solutions assuming long-term (i.e., T min = 100 years) and medium-term (T min = 60+ years) protection of caribou habitat. To estimate the total amount of connected habitat in the absence of harvest, we solved the connectivity objective (6) without harvest and then used this estimate to calculate the proportion of connected area in the solutions with harvest using objective Equations (29)-(31).

| General connectivity patterns
We compared the optimal solutions for two distinct sets of scenarios that prioritized harvest and habitat protection. The original problem included 5,322 forest sites, 25.3 K continuous and 208 K binary variables, but the size of the MIP model was reduced significantly after presolve. The no DCHS problem included 21.3 K continuous and 78.7 K binary variables after presolve and reached an optimality gap below 0.5% in less than 6 hr. The static DCHS problem included 21.3 K continuous and 47.4 K binary variables after presolve and reached an optimality gap of 0.5% within 18-24 hr. The dynamic DCHS problem included 21.3 K continuous and 80 K binary variables after presolve and reached an optimality gap of 0.12%-0.25% after 72 hr. The maximum level of sustainable harvest with given assumptions and constraints approached 0.31 M m 3 -yr −1 in no DCHS scenarios. At this level, suitable habitat can be connected in a long-term network that covers approximately 56% of the Whitesand-Armstrong FMU area. Both harvest priority and habitat protection priority solutions show the bulk of the harvest allocated in close proximity to a network of logging roads in the south-central portion of the FMU (Figure 5a,b). Prioritizing habitat connectivity increases the area of connected habitat by 5.0%-9.5% (Table 4). The habitat protection solutions show unharvested corridors connecting habitat areas on the northern shores of Lake Nipigon and habitats in the undeveloped northern parts of the FMU (Figure 5b). These solutions also avoid allocating harvest in areas with the most suitable category 1 habitat and seasonal use category 2 habitat, instead allocating much of the harvest to eastern parts of the FMU ( Figure 5). Although the harvest allocations in the habitat protection solutions keep more habitat connected than the harvest priority solutions, they increase the mill gate timber cost per unit by as much as $5.4 m −3 (Figure 6). The average timber cost increase was $3.3 m −3 , which is comparable with the size of royalties paid by forest companies for harvesting timber on Crown lands in Ontario. Figure 6b also indicates that the impact of caribou habitat protection on timber supply costs is appreciable even at relatively low harvest levels (0.1 M m 3 -yr. −1 and above). This is because the sites with the most productive mature coniferous stands represent both a desirable source of timber and highly preferred caribou habitat. Their protection moves harvest to lower-productivity sites that are more geographically remote, which increases the timber unit cost.
Imposing habitat connectivity constraints increases the proportion of the study area where forests are harvested a single time over the planning horizon (Table 4). This stems from the allocation of harvest to less productive locations to protect prime caribou habitat. When  connectivity constraints are in place, increasing the harvest volume target requires harvesting larger numbers of these low-productivity sites in addition to the intensively managed area in the south-central part of the FMU, nearly all of which is harvested twice over the planning horizon, regardless of scenario. Overall, an efficient habitat protection strategy is to protect areas with high-use caribou habitat near the northern, southwestern and eastern FMU borders while increasing harvest intensity in the south-central region, which has the lowest hauling costs.

| Impact of DCHS polices
Imposing DCHS rules concentrates harvest within a set of regions that guarantees a sufficient area of undisturbed habitat around these regions. The impact of DCHS is most noticeable in harvest priority scenarios, in particular in the dynamic DCHS scenario (Figure 5e), which imposes a strict penalty on harvesting in neighbouring DCHS regions during the same period. The impact of DCHS is less distinguishable in habitat priority solutions (Figures 5d and 5f) because the creation of large contiguous regions of connected habitat must be offset by allocating harvest to much of the remaining FMU area, regardless of the DCHS rules. With respect to the harvest priority scenarios, applying the static DCHS rules leads to only a small increase of the timber costs compared to no DCHS solutions ( Figure 6 and Table 4). This F I G U R E 6 Impact of the timber harvest target on the area of connected habitat and the timber cost: (a) connected proportion of the total range area, %, vs. timber volume harvest target, million m 3 -year -1 ; (b) mill gate timber price, $ m −3 , vs. timber volume harvest target, million m 3 -year −1 . Solid lines depict harvest priority scenarios and dotted/dashed lines depict habitat connectivity priority scenarios. DCHS, Dynamic Caribou Harvesting Scheduling is because a significant portion of timber harvest occurs in the south-central part of the FMU, which is exempt from the DCHS rules ( Figure 3, callout I) and therefore provides flexibility when allocating harvest. However, in the dynamic DCHS scenarios, which feature a strict penalty on harvesting in adjacent DCHS regions, the cost of timber supply increases sharply and approaches the timber costs in the habitat protection solutions. Still, long-term habitat protection (i.e., T min = 100) imposes a higher premium on the timber supply cost than in harvest priority scenarios with dynamic DCHS because it requires more substantial spatial reallocations of harvest to protect a sufficient amount of prime caribou habitat (see the connected habitat patterns in Figure 5e vs. Figure 5b,d,f). Primarily, the impact of DCHS in habitat priority solutions is captured as changes in the harvest patterns. Harvest tends to be concentrated in DCHS regions proximal to the network of access roads, similarly to the scenario without DCHS, so the coarse-scale harvest pattern, including the total area harvested, remains relatively stable. At low harvest levels, in habitat priority solutions, the cost of timber in the dynamic DCHS scenario is lower than in the static and no DCHS scenarios (Table 4). As mentioned earlier, this is because the stringent constraints of the dynamic DCHS scenario on harvesting in DCHS regions force the model to locate as much harvest as possible in the exempt south-central region (Figure 3, callout I). Doing so allows the most flexible harvesting and decreases the timber cost. While this behavior persists at high harvest levels, the larger volume of harvest cannot be allocated completely to the south-central region, so a sizeable share of the harvest must be apportioned among other DCHS regions. The net result is that, at the highest harvest target (0.3 M m 3 -year −1 ), the timber costs converge at just over $35 m −3 for the habitat priority scenarios with static, dynamic or no DCHS rules. In fact, the costs converge between these scenarios and the harvest priority scenario with dynamic DCHS.
Our results also confirm that the DCHS rules are effective in creating undisturbed habitat space around the harvested regions. For example, in static DCHS solutions the number of adjacent DCHS regions harvested in the same time period is reduced significantly compared to no DCHS solutions, and is further reduced to minimum adjacency levels in the dynamic DCHS solutions (Figure 7a). Habitat priority solutions impose stricter spatial constraints on harvest and usually have slightly more harvest in adjacent regions than the harvest priority solutions (Figure 7a).
In our optimal solutions, the protection of high-use habitat is mandatory but the protection of seasonal habitat is optional. Figure 7b indicates that the habitat priority solutions protect significantly more areas of seasonal habitat than the harvest priority solutions. With respect to the DCHS scenarios, a notable reduction in the area of protected seasonal habitat only occurred in the dynamic DCHS scenarios (Figure 7b). The dynamic DCHS rules have a big impact when the harvest volume target is high and approaches the maximum sustainable harvest limit, otherwise the choice of the objective (i.e., harvest vs. habitat protection priority) has a bigger impact on the area of the protected seasonal habitat than the implementation of DCHS rules.

| Long-term versus medium-term habitat protection
We compared long-term protection solutions, where the minimum habitat protection period T min was 100 years, with solutions emphasizing medium-term protection, with T min set to 60 years or longer. Relaxing the long-term habitat protection requirement allows harvest to extend over a larger area. In general, medium-term habitat protection solutions increased the area of forest stands harvested only once (Figure 8 and Table 5). The total area harvested over the planning horizon T increased, but the hotspots where stands were harvested twice remained the same (Figure 9). Scenarios with the medium-term habitat protection objective connected less habitat area over the long term than scenarios with the long-term habitat protection objective (Table 5). Overall, the solutions with the long-term habitat protection objective (T min = 100 years) enabled protection of, on average, 11.7% more habitat area over the long term while giving up the medium-term (60 year) protection of only 5.4% habitat area (Table 5).
Relatively small geographical differences between the solutions with the long-term and medium-term habitat protection targets can be attributed to constraint (12), which maintains the minimum uninterrupted habitat protection period T min . Possible habitat protection options to maintain a 60-year protection period are limited to deferral of harvest to 40 years or immediate harvest at t = 0 in sites which can regain suitable habitat status in 40 years. This significantly reduces the potential number of sites where such measures could be implemented, and so does not change the location of major habitat protection hotspots (Figure 9).

| DISCUSSION
Incorporating the concept of long-term habitat connectivity into forest planning can help reduce the negative impacts of harvest activities on caribou populations. In our study, the main  regions where long-term protection of habitat was cost-effective included prime habitat areas near the borders of the FMU, including along the north shore of Lake Nipigon. Habitat protection in these regions can be achieved by relocating harvest to the south-central region with a dense network of access roads, and by using more intensive harvest regimes over the planning period. Prioritizing habitat connectivity leads to a small increase in the overall harvest area because harvest has to be allocated to less productive sites to protect prime caribou habitat containing old conifer stands.
Our results indicate that the implementation of DCHS, which is the current harvesting policy in Ontario's boreal forests when caribou populations are present, causes only a minor increase of the timber supply cost in our study area. The low impact of DCHS on the timber price is a result of the decision to exclude the south-central portion of the FMU from the DCHS rules (Figure 3, callout I). Notably, several other FMUs in northern Ontario (such as Black Spruce Forest, English River Forest, Kenogami Forest, and Whiskey Jack Forest) are similar in that respect and have DCHS regions delineated in the northern parts, while southern portions are exempt from DCHS harvesting. The impact of DCHS on the timber supply cost would be more significant if the entire FMU area was subject to DCHS. Strict enforcement of DCHS with the control of harvesting in adjacent DCHS regions increases the cost of timber, but the cost increase does not exceed the cost premiums in the solutions prioritizing protection of caribou habitat. In our study area, the complex configuration of the access road network amidst caribou habitat areas, some of which are high-use, translates to strict harvest limitations when prioritizing the protection  of habitat. The impact of these spatial restrictions on harvest patterns is greater than the impact of DCHS rules currently planned for the FMU. This translates to a higher cost of timber supply, but enables protecting, on average, 7.4% more of the total range area and 12.1% more area with seasonal habitat than in the harvest priority scenarios. In our case, adopting a DCHS strategy in conjunction with the habitat protection objective achieves a (d) (f) F I G U R E 9 Examples of harvest and habitat connectivity patterns for the scenarios with the long-term (T min = 100 years) and medium-term (T min = 60+ years) habitat protection objectives. The maps show the habitat protection priority solutions with the harvest volume target 0.3 million m 3 -yr −1 . No-DCHS scenarios: (a) Long-term habitat protection (T min = 100 years); (b) Medium-term habitat protection (T min = 60+ years). Static DCHS scenarios: (c) Long-term habitat protection (T min = 100 years); (d) Medium-term habitat protection (T min = 60+ years). Dynamic DCHS scenarios: (e) Long-term habitat protection (T min = 100 years); (f) Medium-term habitat protection (T min = 100 years). Dark-shaded areas indicate the habitats connected over 100 years, medium-shaded areas-the habitats connected over 60-less than 100 years and light-shaded areas indicate the habitats connected over less than 60 years. Small and large dot symbols show the forest sites with one and two consecutive harvests over the planning horizon T. DCHS, Dynamic Caribou Harvesting Scheduling reasonable balance between the goals of meeting a desired harvest target and long-term habitat protection, despite the increase in timber supply cost.
Notably, the proposed model was not designed to maintain connectivity corridors between core habitat areas. While the optimal solutions established some corridors between areas with high-use habitat, this was unintentional. Potentially, a set of corridor connectivity constraints could be added to maintain corridor space between high-use habitat areas in different parts of the range, like in the models of St. John et al. (2016) and Conrad et al. (2012). Such a formulation could be useful for planning harvest in areas with isolated caribou ranges, which would require maintaining corridors to facilitate the movement of the animals between the ranges. Habitat connectivity criteria and corridor design are likely to become more important in the future, as the total amount of intact forest suitable to support caribou populations in the Canadian boreal region is likely to decline (EC, 2011).

| Potential model extensions
The current model formulation did not track costs of constructing logging roads, but instead factored such costs into the total tree-to-mill delivery cost value from a particular forest site. Potentially, the location and timing of road construction could be incorporated via a network flow model similar to the formulation presented in Yoshimoto and Asante (2018). This may increase the cost of moving the harvest to remote DCHS regions, which currently have no access roads.
Our model did not impose limits on the minimum width of the corridors connecting the sites with large amounts of habitat. While the optimal solutions delineated sufficiently large contiguous clusters of suitable habitat, some corridors connecting these large areas had a width of a single map cell. In our study, the minimum width of the habitat corridor was set by our chosen spatial data resolution. The choice of a 1-km resolution followed caribou recovery guidelines from Environment Canada (2011) which calls for a minimum 500-m buffer between protected sites and human disturbances (so a point surrounded by a 500-m buffer suggests a minimum 1-km spatial resolution to model habitat connectivity). Commonly, harvest scheduling is performed at finer spatial scales, but for the sake of practicality, it may be necessary, as in our case, to coarsen the resolution of harvest planning while maintaining the minimum width of the established corridors. Potentially, an approach similar to that presented in St. John et al. (2018) could ensure the minimum width of habitat corridors, but is likely to increase the size of the optimization problem substantially.

ACKNOWLEDGMENTS
The funding for this study was provided by Natural Resource Canada's Cumulative Effects Program. Our sincere thanks to Josie Hughes for advice on the caribou habitat model and DCHS and John Pedlar for useful comments on an early version of the manuscript.