Advancing contiguous environmental land allocation analysis, planning and modeling

ABSTRACT

The identification or delineation of a spatially contiguous area within a larger region is an important land use planning problem and is utilized across a range of disciplines and management contexts. This involves the selection of individual spatial units, such as parcels, forest stands, management areas, etc., for an intended goal, purpose, usage and/or activity with constraints on maximum size along with contiguity requirements. The goal of this paper is to present a new spatial optimization model that can be used to simultaneously delineate single or multiple patches, where each is a contiguous set of individual spatial units, meets specific size requirements and is spatially distinct. We structure this new model as a mixed integer-linear program. We also present an application of this new model to a land management problem encountered by the USDA Forest Service involving wildfire mitigation efforts.

KEYWORDS:

• Land acquisition
• patch delineation
• spatial optimization
• wildfire
• contiguity

Introduction

A fundamental land use planning problem involves the selection of spatial units, such as land parcels, forest stands, management areas, and the like, for an intended goal, purpose, usage and/or activity. If each spatial unit is large enough for the intended purpose, then the problem can be quite simple, involving the evaluation of individual units and picking that unit that best meets the desired size, fits within a budget, and other desired criteria. But when all or most of the spatial units are relatively smaller in size, then the problem is substantially more complicated, requiring an exhaustive search for the group of spatial units that form a contiguous area or ‘patch’ to best meet desired objectives and constraints (Shirabe, 2005; Wright et al., 1983). Applications have involved geographical site selection (Cova & Church, 2000; Xiao et al 2002.,), land acquisition (Diamond & Wright, 1991), biological reserve design (McDonnell et al., 2002; Wang & Önal, 2016), viable habitat patch identification for species protection (Church et al., 2003), forest management to reduce potential fire severity and size (Church et al., 2015; Kreitler et al., 2019), and many others (see, Brookes, 2001; Eastman et al., 1995; Malczewski, 2006; Minor & Jacobs, 1994). While there are many more examples that could be offered, of importance here is recognizing the common characteristics shared by many of these problems, which is the identification of the best area(s) for the associated activities, whether it be waste processing, species protection, wildfire fuel treatment or any other context. In a formal sense, one seeks the best collection of spatial units forming a sufficiently large area for its intended use. Additionally, given that the usage or activity is integrated and coordinated, there is often a need for the selected spatial units to be contiguous, or spatially connected (Williams et al., 2004). For example, to support a mated pair of the San Joaquin Kit Fox in California, a contiguous patch of suitable habitat of at least 390–520 ha (1.5 to 2 square miles) is needed (Church et al., 2003). Habitat protection of the Kit Fox therefore requires identifying and selecting suitable patches that best satisfy the needs of individual mated pairs.

The objective of this paper is to present a new model that can be used in various problem settings where the goal is to delineate a prespecified number of project areas, referred to commonly as patches. Each patch needs to be a contiguous set of spatial units, meet size requirements, and be composed of land units that are highly suitable for a particular management action (e.g. species protection, vegetation treatment, etc.). Although this new model is quite general and can be applied in different problem settings, it has been developed and applied to identify the best patches for targeting risk mitigation efforts associated with wildfire.

The risk of wildfire has increased dramatically due to development, fire suppression, selective logging of large fire-resistant trees, and climate change. While fire is a natural phenomenon and historically an important component of the ecosystems in which it occurs, many factors have contributed to increased risks and more serious outbreaks of wildfire – many due to shaping of the environment by humans. Among these factors are changing landscape management practices, including a combination of forest harvesting, land clearing, and fire suppression. The result is the removal of drought tolerant and fire resilient vegetation as well as proliferation of overgrown vegetation. Beyond this, climate change and prolonged drought have increased stress on and mortality of trees and vegetation as well as made them more susceptible to disease and insect infestations. The result has been a build-up of fuels that greatly increases wildfire risk and severity. However, strategic planning of treatment activities, like prescribed burning, harvesting, tree thinning, and mechanically removing ladder fuels and litter, can reduce wildfire potential and shift landscape structure and composition to within the natural range of variation that has proven to be more resilient to wildfire. Indeed, such treatment actions are seen as key to sustainability efforts, recognizing that wildfire is part of a coupled natural-human system in which humans are not only impacted by but also shape the physical environment (Broitman, 2020; Moritz et al., 2014; Munroe et al., 2019; Spies et al., 2014). Wildfires as a phenomenon are closely intertwined with and influenced by the particular ways in which humans have reformed the environment in which they occur and it is clear that neither humans nor wildfires are leaving this coupled system. Ultimately, there is increasing need to embrace co-existence with wildfire (Moritz et al., 2014) and to direct resources to activities that enhance wildfire resilience (Kreitler et al., 2019).

It is precisely this wildfire fuels treatment context that is of interest in this paper, where the best patches to undertake risk reduction measures are sought. Operational considerations employed by the USDA Forest Service, as an example, suggest that it is preferable to have a treatment area of no larger than 40 ha. (100 acres) and contiguous in spatial extent in order to support management and implementation efforts, facilitating the use of needed equipment and enabling personnel access. The goal, therefore, is to delineate contiguous patches within a region that have the greatest potential for reduction of wildfire risk and vulnerability. Finney (2001) has demonstrated that strategic treatment of about 25% of the landscape will significantly enhance resiliency across a landscape. This means that efforts to reduce fuels within selected patches can substantially reduce risk of wildfire across an entire region.

Modelling land use allocation

In a broad sense, the patch delineation problem can be thought of as a form of land use planning, where specific parcels, polygons or land management units (LMUs) are selected for an activity and others are not. Land use planning is fundamentally important (An et al., 2021; Karakostas, 2016). There are many different land use planning models, invariably structured to address individual needs and concerns (Hopkins et al., 1978; Karakostas, 2016; Ligmann‐Zielinska et al., 2008; Minor & Jacobs, 1994; Xiao & Murray, 2019; Yao et al., 2019), and they help to shape the theoretical foundations of land use sought by Turner et al. (2020) and others. Land use planning is inherently concerned with the shaping of a coupled natural-human system and many tools consider this explicitly, including models of forest vegetation cover (Munroe et al., 2019). The regional scale, which includes the forest-level landscape, is of particular relevance to coupled natural-human systems because it includes many types of human and environmental interactions (Wu, 2019). Land use planning models are critical for efforts to improve coupled natural-human systems.

In the context of land allocation, there are two fundamental model structures that underlie most approaches (Church & Murray, 2009): threshold and knapsack. The threshold structure generally seeks to minimize costs subject to the selection of a specified number of units (or equivalently a minimum size). Examples include the approaches of Wright et al. (1983), Gilbert et al. (1985), Diamond and Wright (1991), and Davis et al. (1996), and Shirabe (2005), but a threshold structure is also implicit in Toth et al. (2009). In contrast to a threshold, the knapsack structure seeks to maximize benefit while limiting the number of units that can be selected Examples include Gilbert et al. (1985), Diamond and Wright (1991), and Davis et al. (1996) as well as the primary approach of Shirabe (2005).

Of particular interest in this paper is the knapsack model structure due to its relevance in the context of wildfire mitigation. Ager et al. (2013, 2021) provide an extensive discussion of the knapsack structure for addressing wildfire concerns, indicating that patches are sought for treatment that adhere to a maximum size in order to keep them manageable. In general terms, the knapsack approach has been the subject of much interest and study, with formal discussion and solution details found in Dantzig (1957), among others. Define the following notation:

$\mathit{i}=$ index of spatial units, designated as LMUs (land management units)

${\beta }_i=$ expected benefit of LMU $\mathit{i}$ being selected as a member of a patch

${\alpha }_i=$ area of LMU $\mathit{i}$

$\mathit{T}=$ patch threshold for maximum treated area

$X_i=\left\{\begin{array}{c}1\mathrm{ifLMU}\mathit{i}\mathrm{isselectedfortreatment}\\ 0\mathrm{otherwise}\end{array}\right.$

With this notation, the formulation of the knapsack model is as follows:

Maximize

(1) $\sum _i{\beta }_i{X}_i$(1)

Subject to

(2) $\sum _i{\alpha }_i{X}_i\le \mathit{T}$(2)

(3) $X_i=\left\{0,1\right\}\mathrm{\forall i}$(3)

The objective (1) seeks to maximize total benefit in the selection of LMUs. Constraint (2) stipulates a maximum combined size of the selected LMUs. Binary restrictions are indicated in constraints (3).

This model, (1)-(3), is essentially the formulation given in Ager et al. (2013, 2021) for patch identification, with a note that selected units must form a contiguous patch. No mathematical details are provided for ensuring contiguity, however, as this is accomplished within a heuristic solution process. This is problematic because exact specification and solution are critical for ensuring that the best outcomes possible are considered.

Imposing contiguity restrictions within a patch

The basic problem faced by Ager et al. (2013) and others is that of imposing contiguity between spatial units selected for wildfire mitigation. A patch, or project area, is comprised of smaller land units, or LMUs. The definition of contiguity is absolute in this case, where a patch is contiguous if it is possible to travel between any two members of the patch by moving across only patch members (Church & Murray, 2009; Girvetz & Greco, 2007). There have been a number of approaches to explicitly impose contiguity in land acquisition. Gilbert et al. (1985) and Diamond and Wright (1991) detail implicit enumeration approaches capable of identifying a contiguous patch. An issue, however, is that enumeration-based approaches are computationally limited in the number of spatial units that can be considered. Brookes (1997) proposes a region growing heuristic to ensure contiguity while maximizing the utility of the identified patch. A seminal development was that of Cova and Church (2000), detailing explicit constraints to impose contiguity within a mathematical optimization model, structured based on network flow between selected units. This was effectively extended by Shirabe (2005) to address patch identification for land allocation. There are other formal model structures that have been proposed to enforce contiguity (e.g. Aerts et al., 2003; Toth et al., 2009; Williams, 2002), but unfortunately they are not easily extended to address multi-patch selection.

In the next section, we present a new land use acquisition and allocation model for identifying the best set of patches that satisfy both size and contiguity requirements. The model is formulated mathematically as a mixed integer-linear program. An application of this model is then detailed in the context of wildfire fuel treatment planning for the USDA Forest Service. The results are presented for a study area in a National Forest in California to highlight capabilities of the developed model and demonstrate insights possible, but also to highlight important future research needs. The paper ends with discussion.

Materials and methods

A spatial optimization model has been developed for multi-patch delineation in order to support land allocation associated with wildfire mitigation efforts. This model is formulated as an integer linear programming problem. We assume here that the landscape is divided into many spatial units, called land management units (LMUs). Each LMU is of a known size with a derived benefit for being selected as a part of a patch in terms of the increase in resilience to wildfire that would result if it is treated. Each patch must be composed of a contiguous set of LMUs and not exceed a stipulated total size. The overall objective is to maximize the total value of the selected LMUs across all the delineated patches. Computational processing was done on a desktop personal computer (Intel Xeon E5 CPU, 2.30 GHz with 96 GB RAM) running Windows. XPRESS (version 8.12) via Workbench was used to implement and solve the model. Further details are as follows.

Optimizing single patch selection

Shirabe (2005) developed a land allocation model designed to select individual spatial units that are contiguous. Consider the following additional notation:

$\mathit{M}=$ large number

$N_i=$ set of units that are adjacent to or neighboring unit $\mathit{i}$

$V_i=\left\{\begin{array}{c}1\mathrm{ifunit}\mathit{i}\mathrm{isselectedaspatchsink}\\ 0\mathrm{otherwise}\end{array}\right.$

$Y_{\mathit{ij}}=$ accumulated flow from unit $\mathit{i}$ to $\mathit{j}$ destined for sink

The introduced decision variables are used to identify a patch sink that serves as a collector of flow from all selected spatial units in the patch, where one unit of flow is created when it is selected. In network theoretic terms, each spatial unit can be represented by a node and arcs reflect adjacent or neighbouring units (e.g. $N_i$). Viewed in this way, the significance of such an approach is that network flow is only possible between selected neighbouring units, thereby ensuring contiguity among patch members. The associated knapsack-based single contiguous patch land allocation problem is as follows:

Maximize

(4) $\sum _i{\beta }_i{X}_i$(4)

Subject to

(5) $\sum _i{\alpha }_i{X}_i\le \mathit{T}$(5)

(6) $\sum _{j\in N_i}{Y}_{\mathit{ij}}-\sum _{j\in N_i}{Y}_{\mathit{ji}}\ge X_i-\mathit{M}{V}_i\mathrm{\forall i}$(6)

(7) $\sum _i{V}_i=1$(7)

(8) $\sum _{j\in N_i}{Y}_{\mathit{ij}}\le \left(\mathit{M}-1\right)X_i\mathrm{\forall i}$(8)

(9) $V_i\le X_i\mathrm{\forall i}$(9)

(10) $\begin{array}{rl}& X_i\in \left\{0,1\right\}\mathrm{\forall i}\\ & V_i\in \left\{0,1\right\}\mathrm{\forall i}\\ & Y_{\mathit{ij}}\ge 0\mathrm{\forall i},\mathit{j}\end{array}$(10)

As previously described, objective (4) seeks to maximize total benefit in the selection of spatial units. Constraint (5) stipulates a maximum size of the set of selected units. Constraints (6) indicate a type of conservation of flow for selected, non-sink units (or nodes). Constraint (7) requires that one selected unit be designated as the sink. Constraints (8) link unit selection and flow, limiting the flow out of a unit unless it is selected as part of the patch. Constraints (9) impose that a sink must be a member of the selected units. Binary and non-negativity restrictions are indicated in constraints (10).

A new multi-patch selection model

The model introduced by Shirabe (2005), (4)-(10), is important and noteworthy in many ways, but is limited to the selection of only one patch. As our intent is to identify multiple independent patches simultaneously, an alternative approach is necessary. In order to develop an extension that can be utilized to identify multiple patches, the following additional notation is introduced:

$\mathit{k}=$ index of patches

$\mathit{p}=$ number of patches to be identified

$X_{\mathit{ik}}=\left\{\begin{array}{c}1\mathrm{ifunit}\mathit{i}\mathrm{isselectedfortreatmentaspartofpatch}\mathit{k}\\ 0\mathrm{otherwise}\end{array}\right.$

$V_{\mathit{ik}}=\left\{\begin{array}{c}1\mathrm{ifunit}\mathit{i}\mathrm{isselectedasasinkinpatch}\mathit{k}\\ 0\mathrm{otherwise}\end{array}\right.$

$Y_{\mathit{ijk}}=$ accumulated flow from unit $\mathit{i}$ to $\mathit{j}$ destined for sink in patch $\mathit{k}$

The purpose of this additional notation and associated decision variables is to distinguish unit membership by patch, which enables tracking of patch constraints and conditions, namely maintaining maximum area limitations as well as contiguity within a patch. Given a stipulated number of patches to be identified, a mixed integer-linear program is introduced, referred to the multi-patch selection model:

Maximize

(11) $\sum _i\sum _{\mathit{k}=1}^p{\beta }_i{X}_{\mathit{ik}}$(11)

Subject to

(12) $\sum _i{\alpha }_i{X}_{\mathit{ik}}\le \mathit{T}\mathrm{\forall k}$(12)

(13) $\sum _{j\in N_i}{Y}_{\mathit{ijk}}-\sum _{j\in N_i}{Y}_{\mathit{jik}}\ge X_{\mathit{ik}}-\mathit{M}{V}_{\mathit{ik}}\mathrm{\forall i},\mathit{k}$(13)

(14) $\sum _i{V}_{\mathit{ik}}=1\mathrm{\forall k}$(14)

(15) $\sum _{j\in N_i}{Y}_{\mathit{ijk}}\le \left(\mathit{M}-1\right)X_{\mathit{ik}}\mathrm{\forall i},\mathit{k}$(15)

(16) $V_{\mathit{ik}}\le X_{\mathit{ik}}\mathrm{\forall i},\mathit{k}$(16)

(17) $\sum _{\mathit{k}=1}^p{X}_{\mathit{ik}}\le 1\mathrm{\forall i}$(17)

(18) $\begin{array}{rl}& X_{\mathit{ik}}\in \left\{0,1\right\}\mathrm{\forall i},\mathit{k}\\ & V_{\mathit{ik}}\in \left\{0,1\right\}\mathrm{\forall i},\mathit{k}\\ & Y_{\mathit{ijk}}\ge 0\mathrm{\forall i},\mathit{j},\mathit{k}\end{array}$(18)

Objective (11) seeks to maximize total benefit in the selection of units, consistent with objectives (1) and (4) detailed previously. Constraint (12) stipulates a maximum size for each set of selected units that form an individual patch. Constraints (13) indicate a type of conservation of flow for selected, non-sink units in each patch. Constraints (14) require that only one selected unit be designated as the sink in each patch. Constraints (15) link unit selection and flow, limiting the flow out of a unit unless it is selected for inclusion in a patch. Constraints (16) impose that a sink in a patch must be a member of the selected units for that patch. Constraints (17) indicate that a unit can be included in at most one patch. Binary and non-negativity restrictions are indicated in constraints (18).

Important model properties

The multi-patch selection model, (11)-(18), is actually quite complex in that it involves an expanded form of a multiple knapsack problem with a contiguity requirement for each knapsack. The contiguity conditions (13)-(17) represent a multi-dimensional set that enforces contiguity within each identified patch. In what follows, constraint performance and validity in relation to the planning problem goals are demonstrated and proven.

In order to illustrate how contiguity is imposed in (13)-(16), as well as in (4)-(10), consider the selected units forming a contiguous patch in Figure 1. Focusing only on the selected units, Figure 2 shows labelled nodes along with directed arcs with flow between selected units going to the sink, constituting a feasible solution. Each selected unit (or node) generates one unit of flow in the shown system. Thus, Figure 2 indicates that unit I is the sink, with accumulated flow coming directly from units F, G and H. Specifically, a flow of six comes from unit F, a flow of one from unit G and a flow of one from unit H. The reason for this is that flow from units A, B, C, D and E goes through unit F by way of unit C. If we examine constraints (13) for this solution, assuming this corresponds to patch $\mathit{k}$ with all decision variables $Y_{\mathit{ijk}}$, $X_{\mathit{ik}}$ and $V_{\mathit{ik}}$ equal to zero unless otherwise noted, the following conditions can be observed for selected units (three representative situations are shown, involving units A, C and I, respectively):

(19) ${\sum _{\mathit{j}\in \mathit{N}}}_{{ }_A}{Y}_{\mathit{Ajk}}-\mathit{\hspace{0.17em}}\sum _{j\in N_A}{Y}_{\mathit{jAk}}\ge X_{\mathit{Ak}}-\mathit{M}{V}_{\mathit{Ak}}=1-0\ge 1-\mathit{M}\times 0\left(Y_{\mathit{ABk}}=1\mathit{and}{X}_{\mathit{Ak}}=1\right)=1\ge 1$(19)

(20) $\sum _{j\in N_C}{Y}_{\mathit{Cjk}}-\sum _{j\in N_C}{Y}_{\mathit{jCk}}\ge X_{\mathit{Ck}}-\mathit{M}\mathit{\hspace{0.17em}}{V}_{\mathit{Ck}}=5-4\ge 1-\mathit{M}\times 0(Y_{\mathit{CFk}}=5,Y_{\mathit{BCk}}=2,Y_{\mathit{DCk}}=2\mathit{and}{X}_{\mathit{Ck}}=1)=1\ge 1$(20)

(21) $\sum _{j\in N_I}{Y}_{\mathit{Ijk}}-\sum _{j\in N_I}{Y}_{\mathit{jIk}}\ge X_{\mathit{Ik}}-\mathit{M}\mathit{\hspace{0.17em}}{V}_{\mathit{Ik}}=0-8\ge 1-\mathit{M}\times 1(Y_{\mathit{FIk}}=6,Y_{\mathit{GIk}}=1,Y_{\mathit{HIk}}=1,X_{\mathit{Ak}}=1\mathit{and}{V}_{\mathit{Ik}}=1=-8\ge -\mathit{M}+1$(21)

Figure 1. Selected contiguous patch in a region.

Figure 2. Network flow representation among selected units forming a contiguous patch.

All of these conditions are satisfied, making it a feasible solution. If one or more units were selected that did not connect to the sink, the associated constraints would not be valid, and the solution would be infeasible. Thus, the model constraints ensure that flow originates from selected units, and must be destined for the sink, either directly or through other selected units. This structure works because only a selected unit can accept incoming flow from neighboring units, as stipulated in constraints (15). Consider the case for unit A: $\sum _{j\in N_A}{Y}_{\mathit{Ajk}}\le \left(\mathit{M}-1\right)X_{\mathit{Ak}}$. Since unit A is selected, then $X_{\mathit{Ak}}=1$ and only $Y_{\mathit{ABk}}=1$ (all other $Y_{\mathit{Ajk}}=0$). This means that this constraint becomes $1\le \left(\mathit{M}-1\right)\times 1$. Otherwise, all outgoing flow would be forced to be zero.

The above discussion was descriptive in nature, attempting to illustrate what is happening in the model through imposed constraints. The properties of the formulation are now proven. This is done by first making the following observations.

1. Implicit spatial unit feasibility. The three models, knapsack (1)-(3), single patch (4)-(10) and multi-patch (11)-(18), implicitly assume that an individual spatial unit satisfies the knapsack constraint, (2), (5) and (12), respectively. That is, ${\alpha }_i\le \mathit{T}$ for each spatial unit $\mathit{i}$. If this is not the case (e.g. ${\alpha }_i>\mathit{T}$ for one or more units $\mathit{i}$), then clearly constraints (2), (5) and (12) in each model prevent that unit from ever being selected. A number of options exist for dealing with this situation. A spatial unit where ${\alpha }_i>\mathit{T}$ could simply be removed from consideration and not included in the model (or kept in the model, recognizing it will not be selected). Other options may be necessary, however, if the spatial unit is in some way feasible. One is that the unit could be subdivided in such a way that the new spatial units satisfy the limit. That is, define new spatial units that partition $\mathit{i}$, where ${\alpha }_i=\sum _l{\hat{\alpha }}_{\mathit{il}}$ and ${\hat{\alpha }}_{\mathit{il}}\le \mathit{T}$ for each subunit $\mathit{l}$. Another option is to modify the model formulations to allow for the selection of unit $\mathit{i}$ even though it technically violates the limit $\mathit{T}$. This is possible, but the details associated with such an approach are beyond the scope of this paper. The implicit spatial unit feasibility observation is useful in the model validity proofs that follow.

2. There is sufficient capacity across a region for the selection of $\mathit{p}$ patches. Sufficient capacity for the selection of $\mathit{p}$ patches suggests that $\sum _i{\alpha }_i\ge \mathit{T}\times \mathit{p}$. If this is not true, then two conditions result. One is that patches would be selected that are significantly less than the threshold $\mathit{T}$, assuming that the number spatial units is greater than $\mathit{p}$. All spatial units would be selected, making this a trivial and uninteresting optimization problem. The second condition is that the number of spatial units is less than $\mathit{p}$. In this case, all units would be selected, but the number of patches would be less than $\mathit{p}$. Again, this is a trivial case that would not need optimization for analysis, but also suggests that $\mathit{p}$ has been inappropriately set.

3. Most spatial units exhibit a positive benefit for selection. In the context of land use acquisition, there is generally a quantifiable and positive benefit, ${\beta }_i>0$, to the selection of most spatial units. Specifically, if the set of positive benefits units, $\mathrm{\Phi }=\left\{\mathit{i}|{\beta }_i>0\right\}$, satisfies $\sum _{\mathit{i}\in \mathrm{\Phi }}{\alpha }_i\ge \mathit{T}\times \mathit{p}$, then there would incentive to select a unit if possible, though contiguity requirements may necessitate the inclusion of units where ${\beta }_i\le 0$.

Given these three observations, there are two important aspects of the multi-patch selection model, (11)-(18), to prove. This includes the identification of $\mathit{p}$ unique patches and the ability to ensure contiguity within each identified patch.

Lemma 1: The number of identified patches in a multi-patch selection model, (11)-(18), feasible solution will be precisely $\mathit{p}$.

Proof: Assume that only $\mathit{p}-1$ patches are identified and that $X_{\mathit{ik}}=0$ for a case where $\mathit{p}$ patches are supposed to be selected. As most ${\beta }_i>0$ for each $\mathit{i}$ (observation 3), then there is an incentive to always select a spatial unit if at all possible given the orientation of objective (11), and similarly for objective (4). Further, observation 1 establishes individual unit feasibility. The selection of fewer than $\mathit{p}$ patches would imply that for any patch $\mathit{k}$, $X_{\mathit{ik}}=1$ gives $\sum _{\mathit{j}\left(\mathit{j}\ne \mathit{i}\right)}{\alpha }_j{X}_{\mathit{jk}}+{\alpha }_i{X}_{\mathit{ik}}>\mathit{T}$, or equivalently that $\sum _i{\alpha }_i\le \mathit{T}\times \left(\mathit{p}-1\right)$. However, this is a contradiction to observation 2, that $\sum _i{\alpha }_i\ge \mathit{T}\times \mathit{p}$. Thus, the number of identified patches will be precisely. Q.E.D.

The remaining issue to be addressed concerns patch contiguity in the multi-patch section model, (11)-(18), as well as the single patch model, (4)-(10). Constraints (13) and (15) are the major flow components in the model, requiring a spanning tree to connect selected units in a patch. However, constraints (14) and (16) are important as well, requiring each patch to have a sink and restricting a spatial unit from being in at most one patch, respectively. A theorem regarding network flow now follows.

Theorem 1: Flow constraints (13)-(16) impose patch contiguity in the multi-patch model, (11)-(18).

Proof: Without loss of generality, consider a given patch $\mathit{k}$. Assume that a particular spatial unit $\mathit{i}$ has been selected, $X_{\mathit{ik}}=1$, but it is not contiguous or connected to other selected units in patch $\mathit{k}$. This implies there is one or more additional units selected as well. Consider first the case of an additional single selection, unit $\mathit{j}$ (e.g. $X_{\mathit{jk}}=1$). If unit $\mathit{i}$ is a sink ($V_{\mathit{ik}}=1)$, then by the corresponding constraint (13) with $X_{\mathit{ik}}=1$, we have $\sum _{j\in N_i}{Y}_{\mathit{ijk}}-\sum _{j\in N_i}{Y}_{\mathit{jik}}\ge 1-\mathit{M}\times 1$. Since unit $\mathit{i}$ is assumed not connected, this implies $\sum _{j\in N_i}{Y}_{\mathit{ijk}}=0$ and $\sum _{j\in N_i}{Y}_{\mathit{jik}}=0$. Thus, constraint (13) becomes $0-0\ge 1-\mathit{M}$ (or $0\ge 1-\mathit{M}$ or $\mathit{M}\ge 1$). The corresponding constraint (15) would be $\sum _{j\in N_i}{Y}_{\mathit{ijk}}\le \left(\mathit{M}-1\right)\times 1$ as $X_{\mathit{ik}}=1$. Since $\sum _{j\in N_i}{Y}_{\mathit{ijk}}=0$ because it is assumed unconnected, then this constraint becomes $0\le \left(\mathit{M}-1\right)$. However, since another unit is assumed selected (e.g. $X_{\mathit{jk}}=1$) and patch $\mathit{k}$ can only have one sink by constraint (14), this means that $V_{\mathit{jk}}=0$. Further, all other potential sinks in patch $\mathit{k}$ are not possible as a pre-condition in constraints (16), so $V_{\hat{j}\mathit{k}}=0$ for all other units $\mathit{j}$ associated with patch $\mathit{k}$. Thus, the associated unit $\mathit{j}$ constraint (13) is $\sum _{\hat{j}\in N_j}{Y}_{j\hat{j}\mathit{k}}-\sum _{\hat{j}\in N_j}{Y}_{\hat{j}\mathit{jk}}\ge 1-\mathit{M}\times 0$. This simplifies to $\sum _{\hat{j}\in N_j}{Y}_{j\hat{j}\mathit{k}}-\sum _{\hat{j}\in N_j}{Y}_{\hat{j}\mathit{jk}}\ge 1$. But since unit $\mathit{j}$ is the only other selected unit, this implies that $\sum _{\hat{j}\in N_j}{Y}_{j\hat{j}\mathit{k}}=0$ and $\sum _{\hat{j}\in N_j}{Y}_{j\hat{j}\mathit{k}}=0$, or $0-0\ge 1$, or simply $0\ge 1$. This is violation of the set of constraints (13) and is not possible, but more importantly is a contradiction of the assumption. The reverse sink situation ($V_{\mathit{ik}}=0$ and $V_{\mathit{jk}}=0$) leads to the same contradiction. The situation involving three or more selected units in a patch is slightly more involved, but it follows from the two-unit case that constraints (13)-(16) prohibit non-contiguous selection of units, including well known cases of cycles, thereby proving the theorem. Q.E.D.

The significance of the observations, lemma and theorem is that they prove that the multi-patch selection model, (11)-(18), will result in $\mathit{p}$ patches being identified and each patch will form a contiguous area. The application results that follow illustrate that this is indeed true.

Case study data

As indicated previously, the particular land use allocation context investigated in this paper involves wildfire fuels treatment carried out by the USDA Forest Service. The case study consists of five different Potential Operational Delineations (sub-regions) across the Stanislaus National Forest in California, shown in Figure 3. The data is part of the Social and Ecological Resilience Across the Landscape (SERAL) project seeking to increase forest resilience by shifting vegetation conditions to within their natural range of variation. Risk assessments were carried out by the SERAL team, with a focus on the difficulty of wildfire suppression. The LMU represents the smallest area considered for treatment decision-making purposes. LMUs are delineated to maintain relative homogeneity in terms of vegetation type and topography. The objective benefit of treatment is a measure of the departure from what would be a fire resilient area, so the greater the value, the less resilient an LMU would be to wildfire. A value of 0 indicates that an LMU is resilient to wildfire.

Figure 3. Case study areas.

The five study areas consist of 57, 130, 315, 608 and 1343 LMUs (Figure 3). In all application instances, the maximum size of a patch is 100 acres, reflecting a manageable project area for carrying out wildfire treatments. Neighbouring units were considered those that share a common point or edge, the so-called queen criterion. The 57 LMU area has unit sizes ranging from 2.00 to 55.15 acres, with an average size of 13.12 acres. The benefits of wildfire mitigation range from 0.0 to 167.75 (average of 16.59). The number of patches considered was 1, 2 and 3, with three patches resulting in approximately 30% of the area selected for treatment. Such a threshold reflects a state where substantial benefit in wildfire risk reduction is offered to the entire planning area (see, Finney, 2001). The 130 LMU area has unit sizes ranging from 2.00 to 101.85 acres, and an average size of 13.48 acres. The benefits range from 0.0 to 46.71 (average of 7.78). The number of patches considered was 1, 3 and 6, with six patches resulting in 30% of the area selected for treatment. The 315 LMU area has unit sizes ranging from 0.22 to 77.84 acres, and an average size of 10.41 acres. The benefits range from 0.0 to 94.70 (average of 12.34). The number of patches considered was 1, 5 and 11, with eleven patches resulting in approximately 30% of the area selected for treatment. The 608 LMU area has unit sizes ranging from 0.22 to 70.28 acres, and an average size of 9.12 acres. The benefits range from 0.0 to 119.73 (average of 8.83). The number of patches considered was 1, 5, 10, 15 and 19, with 19 patches resulting in approximately 30% of the area selected for treatment. The 1343 LMU area has unit sizes ranging from 0.89 to 217.06 acres, and an average size of 11.75 acres. The benefits range from 0.0 to 198.80 (average of 13.69). The number of patches considered was 1, 5, 10, 15, 20, 30, 40, 50 and 55, with 55 patches resulting in approximately 30% of the area selected for treatment.

Results

The multi-patch selection model was utilized for management planning of each of the five study areas for a range of scenarios. A maximum processing time of 60 minutes was imposed when solving any problem instance. The total benefit in the improvement of wildfire resilience, objective (11), is reported along with the resulting optimality gap (%), if any, as well as total solution time for each of the 23 planning problem application instances. A summary is given in Table 1. Information in Table 1 includes the planning application (referenced by the number of LMUs) and number of patches to be identified. Also provided in Table 1 is the associated number of constraints and decision variables, roughly suggestive of computational difficulty that can be expected, with larger problem instances generally more difficult to solve.

Table 1. Computational results summary using XPRESS.

Application	# Constraints	# Decision variables	Objective, (11)	Optimality gap (%)	Time (sec)
57 (p = 1)	230	410	312.89	0.00	0.5
57 (p = 2)	403	820	506.21	0.00	11.2
57 (p = 3)	576	1,230	667.02	0.00	119.5
130 (p = 1)	522	974	131.24	0.00	185.2
130 (p = 3)	1,306	2,922	347.92	5.19	3600.0
130 (p = 6)	2,480	5,844	575.75	8.16	3600.0
315 (p = 1)	1,262	2,392	206.50	0.01	99.5
315 (p = 5)	5,050	11,960	959.94	7.72	3600.0
315 (p = 11)	10,732	26,312	1903.59	9.62	3600.0
608 (p = 1)	2,434	4,776	211.82	0.00	553.4
608 (p = 5)	9,738	23,880	979.64	5.67	3600.0
608 (p = 10)	18,868	47,760	1776.44	11.30	3600.0
608 (p = 15)	27,988	71,640	2555.36	10.24	3600.0
608 (p = 19)	35,302	90,744	2975.32	14.68	3600.0
1343 (p = 1)	5,374	10,590	284.62	0.00	4.2
1343 (p = 5)	21,498	52,950	1132.167	16.81	3600.0
1343 (p = 10)	41,653	105,900	2115.35	17.30	3600.0
1343 (p = 15)	61,808	158,850	3103.10	14.14	3600.0
1343 (p = 20)	81,963	211,800	3869.37	17.74	3600.0
1343 (p = 30)	122,273	317,700	5435.46	19.31	3600.0
1343 (p = 40)	162,583	423,600	6523.70	26.69	3600.0
1343 (p = 50)	202,893	529,500	6024.03	64.58	3600.0
1343 (p = 55)	223,048	582,450	7501.99	42.62	3600.0

Most problem instances resulted in optimality gaps, indicating that an optimal solution cannot be confirmed within the imposed time limit. Reasons for this are associated with exact solution techniques, in this case mixed integer-linear programming relying on an approach called branch and bound. An optimal solution may not be found due to time limits or that it simply cannot found and/or verified. Nevertheless, solution quality can be established given valid bounds on the solution found, and the optimality gap percent suggests that the reported solution is at least within the indicated theoretical bound.

The 57 LMU application was solved for 1, 2 and 3 patches. Optimal results are reported in each case. Note that problem size grows (number of constraints and number of decision variables) as the number of patches to be identified increases, and as a result solution time (and therefore difficult) also increases. The optimal single patch solution associated with the 57-unit Potential Operational Delineation is shown in Figure 4, illustrating contiguity among selected units that form a patch. The size of this patch is 99.41 acres, maintaining the 100-acre maximum for a patch, giving a total benefit of 312.89. Treating all 57 units would give a total benefit of 945.87, but this is not possible due to budgetary considerations. Thus, the optimal single patch improves resiliency of the Potential Operational Delineation by 33.08%, representing the best possible outcome for improving wildfire resiliency of this region through patch-based mitigation efforts. Of course, this can be increased if the number of patches is increased. Table 1 indicates resiliency benefits of 53.52% (2 patches) and 70.51% (3 patches), respectively, for this region under alternative scenarios.

Figure 4. Selected patch (p = 1) for the 57 LMU application.

The 130 LMU application summarized in Table 1 indicates that as the number of patches is increased from 1 to 6, problem size increases substantially. In this application, optimal results could only be verified for the single patch case. When the number of patches is 3 (p = 3), the total benefit of selected units is 347.92, but a gap of 5.19% remains in the solution process. Thus, this solution can only be characterized as being within 5.19% of the optimum. This highlights that the associated optimization model is challenging to solve. The 3 patch solution is shown in Figure 5. Again, contiguity among patch members can be observed as can the three distinct patches. As stipulated, each patch maintains the 100-acre maximum. Further, increased wildfire resiliency is possible through the treatment of more patches, with total benefits of 131.24 (p = 1), 347.92 (p = 2) and 575.75 (p = 3) indicated in Table 1 for the case of 130 LMUs.

Figure 5. Selected patches (p = 3) for the 130 LMU application.

Similar trends in solution capabilities can be observed for the other planning problem applications, with the 315, 608 and 1343 LMU application results detailed in Table 1 showing that a single patch (p = 1) can be solved optimally, but applications involving more than one patch are substantially more difficult to solve. As a result, optimality gaps remain for all cases except p = 1. Not particularly surprising is that for the larger planning problem instances, larger optimality gaps result. For example, the largest optimality gap reported in Table 1 is 64.58% for the 1343 LMU p = 50 case. Nevertheless, each solution identifies p-patches, and each patch maintains the 100-acre maximum and is contiguous among patch members. Representative solutions for each application include Figure 6 depicting p = 5 patches for 315 LMUs, Figure 7 giving p = 10 patches for 608 LMUs and Figure 8 highlighting the p = 15 patch solution for 1343 LMUs.

Figure 6. Selected patches (p = 5) for the 315 LMU application.

Figure 7. Selected patches (p = 10) for the 608 LMU application.

Figure 8. Selected patches (p = 15) for the 1343 LMU application.

[Figures 6, 7 & 8 near here]

Discussion

This paper discussed land use acquisition and allocation modelling for identifying the best set of patches that satisfy size and contiguity requirements. The multi-patch selection model was introduced and formulated mathematically as a mixed integer-linear program, developed to support planning needs based on important spatial and aspatial characteristics, including a maximum size of areas and the need for selected units to be contiguous. Theoretical properties and characteristics of the model were discussed and proven associated with identifying a pre-specified number of patches and how contiguity is imposed. A wildfire fuels treatment context was described in support of USDA Forest Service management goals and objectives. Application results were presented using five planning areas, highlighting the capabilities of the developed model and demonstrating insights possible. The model developed in this paper is the first to address optimal multi-patch selection. Results from this research have shown that planning efforts can be supported by modelling, enabling important insights to be identified, discussed and considered. Science-based planning efforts are a prerequisite for achieving sustainable coupled natural-human systems, particularly in the context of wildfire risk and vulnerability.

A number of additional discussion items are worth mentioning. One is that the theoretical properties of the multi-patch selection model reflected in Lemma 1 and Theorem 1 are evident in reported results. In particular, the specified number of patches results in that many projects areas being identified, each not exceeding the permissible maximum. Further, each patch is contiguous, consisting of other selected LMUs. The fifteen patches shown in Figure 8 provide an illustration of this, as do Figures 4, 5 and 6.

A second discussion point is that the spatial patterns of selected patches can vary significantly. They many be highly clustered, as in the case of the ten patches shown in Figure 7, or more dispersed, as is the case in Figure 8. An interesting question for coupled natural-human systems, and sustainability more generally, is whether there may be inherent benefits to spatial concentration or dispersion of patch-based wildfire mitigation efforts. Evidence suggests that certain specifies thrive under dispersed activities, while others prefer concentrated activities. Wildfire resiliency and relationships to historic natural conditions make further understanding of patterning of activities important in various ways, with approaches like the multi-patch selection model enabling this to be more systematically studied.

A third discussion item is that high quality, feasible solutions can be found within established quality bounds. Use and application of the multi-patch selection model enables the best wildfire resiliency outcomes to be identified for a pre-specified number of patches, highlighting differences is different treatment possibilities. Summarized in Table 1 was different patch options up to approximately 30% of a region. This is critical in supporting analysis, planning and management efforts. Further, it is computationally viable to be addressed such problems using a commercial optimization solver, particularly in the case of a single patch. It was demonstrated that the developed model can be solved to optimality within a fairly small amount of computational time when identifying the best single patch. There are clear sustainability implications associated with the selection of a patch that offers the greatest wildfire resiliency, making the use of heuristic approaches potentially problematic if they are not of equally high quality. Though not reviewed in the paper, the ForSys heuristic detailed in Ager et al. (2013, 2021) to identify patches are up to 6.69% less efficient (4.56% average across the 23 planning problem instances examined in Table 1) than the exact approach based on the multi-patch selection model. Again, the resiliency and sustainability impacts of such differences are significant.

A final discussion point is that Table 1 highlights that as the number of patches to be selected increases, the difficulty of solving individual problem instances increases significantly. In fact, a majority of the problem instances were not solved to proven optimality within an hour of computer time. Unfortunately, there are some rather significant optimality gaps that remain after one hour of computing. This means that there is room for improvement, either in the way a solver might be employed or in the development of an enhanced model structure. We did not test-specific solver settings or solution strategies, as we relied on default settings. This is one area that needs to be explored. Another is the possibility of defining a tighter form of the contiguity constraint structure or employing a specialized iterative process in solving the multi-patch model (e.g. Duque et al., 2018 for the p-regions problem). Alternative constraints that enforce contiguity might also be tested within the multi-patch selection model (e.g. the spanning tree structure of Duque et al., 2011). In developing the multi-patch selection model, we considered employing the primal-dual graph of Williams (2002) to impose contiguity as computational tests of the primal-dual graph structure have proven to be quite promising. This approach is based upon an assumption that conditions for adjacency/neighbour relationships can be represented as a planar graph. Unfortunately, the adjacency conditions associated with touching at both points and boundary edges (as found in forestry and many other planning contexts) cannot be structured as a planar graph. Consequently, the primal-dual graph approach was not an option. Altogether the results presented here underscore the difficulty in generating optimal solutions for sizable problems. It is worthy to note that the size of the problems solved here are considerably larger than many of the contiguity models considered in previous research (e.g. see, Duque et al., 2011; Shirabe, 2005, 2009). It is also important to underscore the fact that the new model presented here can be used in a number of other problem settings, including land acquisition and biological reserve design.

Acknowledgments

We wish to acknowledge staff of the Stanislaus National Forest, the Social and Ecological Resilience Across the Landscape (SERAL) project team, and Dr. Patricia Manley of the USFS Pacific Southwest Research Station.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Funding for this research through the Pacific Southwest Research Station of the U.S. Forest Service, Department of Agriculture, is gratefully acknowledged (#20-CR-11272139-061);USDA Forest Service Grant;