Assessing the Potential to Increase Landscape Complexity in Canadian Prairie Croplands: A Multi-Scale Analysis of Land Use Pattern

Increasing natural vegetation in agricultural landscapes can create habitat for beneficial organisms such as pollinators and the natural enemies of crop pests. Adding perennial vegetation can also support biodiversity conservation and climate change mitigation objectives. However, implementing such changes to agricultural land use across large geographic areas will require a strategic approach. This study examined the amount and distribution of uncultivated areas in Canadian prairie croplands, focusing on Alberta's agricultural zone (226,543 km2). The aim was to identify locations in this region that have potential for increasing non-crop land cover within fields. This assessment was based on a multi-scale model of landscape complexity that described the distribution of land cover as a function of the distance from field centers. It is based on the assumption that the land cover in the field neighborhood is an informative index of how much non-crop area might realistically be maintained or restored in the field itself; i.e., because neighboring lands will reflect the local environmental conditions that support the growth and establishment of non-crop vegetation as well as the likelihood that crop growers will remove areas from production. The model identified variation across the region in land cover distribution, with regions at latitudes between 52°N and 55°N demonstrating the greatest contrasts in the amount of non-crop land between the field and the field neighborhood scale. These findings suggest that there remains capacity for land use decision-makers to optimize the distribution of non-crop land covers in ways the reduce the differences between these scales (i.e., to increase non-crop covers within fields to better represent the neighborhood proportions). Modeling also revealed scale-dependent patterns, such as field margins without crops (400–500 m from field centers) broadly distributed across the region, and evidence that gradients in moisture and temperature have interacted with land use decisions to shape the proximity of non-crop area to fields.


INTRODUCTION
Increasing agricultural production to meet a growing global demand will require expansion and intensification of croplands (Godfray et al., 2010;Laurance et al., 2014;Egli et al., 2018), but this presents challenges for mitigating climate change and conserving biodiversity (Kleijn et al., 2009;Bustamante et al., 2014;Dalu et al., 2017;Karp et al., 2018). Ecological intensification of agriculture may offer a compromise by leveraging ecosystem services provided by organisms that boost yields while also minimizing impacts on natural systems (Bommarco et al., 2013).
One proposal is to promote non-crop vegetation within or near fields, which would reduce the simplification of agricultural landscapes caused by clearing and field expansion (Landis, 2017). Heterogeneity in non-crop land covers (hereafter landscape complexity), creates a greater number of off-field habitats providing opportunities for "spillover" into the crop (Birkhofer et al., 2018). The beneficial organisms that use these habitats may provide regulating services such as pest control and pollination to the surrounding crop. For example, complex agricultural landscapes have been associated with higher abundance and diversity of birds (Boesing et al., 2017), bats (Ancillotto et al., 2017), flower-visiting insects (Duarte et al., 2018) and the natural enemies of crop pests (Chaplin-Kramer et al., 2011), among other organisms.
There are further arguments for increasing landscape complexity. Restoring or augmenting semi-natural and other non-crop vegetation is an opportunity for carbon storage when trees, shrubs, and other perennial plants are permitted to grow (Smith, 2014;Lamb et al., 2016;Hungate et al., 2017;Williams et al., 2018). It may also support several other regulating ecosystem services through an increase in plant diversity (Asbjornsen et al., 2014), and conserve habitat for organisms that may not have direct benefits to agriculture (Phalan et al., 2011;Tscharntke et al., 2012).
However, meaningfully achieving these "win-win-wins" for agriculture, biodiversity and climate change mitigation will require that land use initiatives be implemented at broad geographic extents. Determining which parts of an agricultural region may be more or less amenable to improving landscape complexity can be used to target early efforts to those areas where there is greatest chance of success. For example, it may be easier to convince growers and land-owners to take land out of production in areas where there is already evidence of higher landscape complexity, but it is not uniformly distributed across all fields. In effect, the landscape context provides a measure of what improvements may be feasibly implemented. This study examine landscape complexity in the Canadian Prairies. The region ranks among the world's largest contiguous agroecosystems, and has a cultivated footprint of more than 0.5 million km 2 (Agriculture and Agri-Food Canada, 2015). The focus is cropland in Alberta, which is distributed across most of the environmental and land use variability of these temperate grasslands (Ecological Stratification Working Group, 1995).
The primary aim is to identify which parts of this region have the highest potential for improving landscape complexity. That is, to identify where introducing additional non-crop land cover into fields and their surroundings may face the fewest challenges to implementation. For example, land conversion decisions may face resistance because natural vegetation is slow or costly to establish in certain regions, or because the land is highly-productive and value is placed on maximizing the area in production.
Finding such areas is done by building a spatial model of the distribution of non-crop land cover and how it varies with proximity to field centers. In this multi-scale approach, distributions of land cover at broad scales provide a target for the introduction of land cover within fields. Secondary aims of this study are to identify patterns in the distribution of these land covers, and how these may be associated with the broad environmental gradients that structure both vegetation and productivity.

Study Area
The geographic focus is Alberta's agricultural zone, one of the most intensively cropped regions in the world (Foley et al., 2005), and an area that represents 30% of Canada's cropland and 20% of its annual farm income (Statistics Canada, 2016; Alberta Agriculture and Forestry, 2017). The dominant land cover is cropland, primarily in a 3-years cereal grain, oilseed, and pulse rotation. Forage lands are a secondary land cover vegetated in both introduced and native perennial grasses that are used for pasture and hay. Other semi-natural areas not in crop production are found both within and adjacent to fields, and they vary in frequency, size and type spatially along environmental gradients. These include patches of perennial grasses, shrubs, and deciduous trees. Planted tree shelterbelts, grass and forb road margins and wetlands of different classes, often surrounded in perennial vegetation, are also common throughout most of the study area (Doherty et al., 2018).

Land Cover
To characterize the variation in land cover across the region, a composite raster map was produced at 30 m resolution by combining data from several published sources in the order listed ( Table 1). The product included a variety of land cover classes which were then simplified thematically to crop and non-crop areas. Non-crop areas also included paved and other humanmodified areas such as roads, buildings, farm yards and in-field oil and gas well pads as these areas are likely to be surrounded in perennial vegetation and therefore contribute to landscape complexity (Forman, 2009). The binary thematic resolution of the land cover map was necessary given that better resolved and accurate products do not exist for this large extent. While this restricts inference about the types of vegetation that may contribute to landscape complexity, it is a simplification that may improve interpretation by reframing land cover as a land use decision (i.e., for crop production or otherwise).

Measuring Landscape Complexity
Landscape complexity was measured from the land cover map using randomly-selected field centers as sampling locations following an algorithmic approach implemented in R. Probable fields were identified using a nineteenth-century land survey which consistently divided the province of Alberta into 259 ha (1 mile by 1 mile) sections (Larmour, 2005). As a consequence of the gridding of the region, fields typically are nested within a section, with the quarter section (square subdivisions; 805 m by 805 m) describing the boundaries of most fields. The field centers ( Figure 1) were randomly selected from all quarter section centroids with the conditions: (1) that the quarter must contain at least 50% crop cover (to ensure it at least partly represents a crop field); and (2) that the centers be no closer than 2.5 km FIGURE 1 | Illustration of the multi-scale method for measuring the proximity of non-crop land covers to field centers. (A) Forty-five annuli of decreasing radii, but approximately equal area, are measured from each field center; (B) Field centers are positioned at the center of quarter sections, the spatial unit that describes most fields in the region, and spaced a minimum of 2.5 km apart; (C) Counts of non-crop area raster cells and the total number of cells in each annulus in A used to find proportional cover. Non-crop raster cells are counted in the smallest annulus in which they fall. apart (to improve statistical independence). These conditions were implemented using an algorithm that iteratively tested a randomly-permuted list of field centers for inclusion until no more could be added.
Measuring landscape complexity requires assessing the variation in landscape structure. Metrics that have been used to describe landscape complexity are typically measured at a certain radius from a focal location, for example, by calculating the area of focal land covers, the proportional composition of crop, or the mean habitat patch size (e.g., Boesing et al., 2017). Because the distance from a field at which a non-crop land cover may have an effect is typically not known, this study uses an alternative approach that creates a multivariate index (or a "curve") integrating the amount of land cover over multiple distances. This landscape complexity curve is determined by finding the proportion of non-crop area in 45 annuli of increasing radii (150-1,006 m) each of which is pre-determined to cover approximately the same total area to minimize sampling bias associated with smaller annuli (Figures 1A,C). Pixels from the rasterized land cover map were recorded in the annulus with the smallest radius in which they fell.
Landscape complexity, as defined here, is therefore a multiscale metric for the proximity of non-crop land covers to the crop and more generally as an index of the variation in non-crop area both within and in the vicinity of a probable crop field. Annuli with the smallest radii describe conditions within the field itself (150-400 m; "within-field scale"), intermediate radii describe the field margins (400-500 m; "field-margin scale") and larger radii which sample from multiple neighboring fields and seminatural areas capture broader landscape variation (500-1,006 m; "neighborhood scale"). The neighborhood scale is intended to capture the local land cover conditions, and is applied in this study as an achievable upper target for restoring or augmenting non-crop land cover within the field itself. Because this scale samples neighboring fields it also measures how much of the landscape is in crop production.

Potential Environmental Drivers of Landscape Complexity
The study region encompasses several environmental gradients that influence both the dominant type of vegetation and primary productivity (e.g., temperature, moisture, latitude, and elevation; Ecological Stratification Working Group, 1995). These gradients also have the potential to affect the frequency of clearing and the rate of natural regrowth of vegetation, and may therefore interact with the behavior of land use decision-makers to shape the observed landscape complexity. Five covariates were chosen to test the importance of these environmental drivers (Figure 2). Mean annual temperature, total annual precipitation and plant available water in 120 cm of soil (April through August) were calculated as 9-years means (2009-2017) from an interpolated data product at the resolution of a survey township (36 sections; 93 km 2 ; Alberta Agriculture and Forestry, 2018). Also included were agricultural limitations for crop production, a classification chiefly based on the workability of the soil by mechanized agricultural equipment (Agriculture and Agri-Food Canada, 2012), and ecoregions, a well-established categorical assessment of contiguous areas with consistent climate characteristics, similar vegetation, soil types, and elevation (Ecological Stratification Working Group, 1995).

Modeling Landscape Complexity at Multiple Scales
Landscape complexity curves were modeled using a type of functional data analysis (FDA) known as function-on-scalar regression (Kokoszka and Reimherr, 2017;Wood, 2017). These models have functions (i.e., the landscape complexity curve) rather than the typical scalar values as response variables. A continuous function Y(x), was estimated from discrete data (e.g., Figure 1C, light line) using the proportion of non-crop area at each annulus radius, x. The model had the general form, where µ (x) is the intercept function (i.e., the mean landscape complexity curve), α p(i) and α p(i) (x) represent the constant effect and functional effect, respectively, of level i of categorical variable p, β q q, x reflects the functional effect of continuous variable q, β r 1 r 2 (r 1 , r 2 , x) is the functional spatial effect given northing, r 1 , and easting, r 2 , and ǫ(x) is the residual function. The model was fit using mgcv, a generalized additive modeling (GAM) package for R (Wood, 2017) using restricted maximum likelihood, by first setting up models with the refund package for R (Goldsmith et al., 2019) which provides optimized settings for FDA. The model assumed a Binomial data-generating process. Categorical constant effects were estimated parametrically, while categorical and continuous functional effects were estimated non-parametrically, with µ(x) and α p(i) (x) functions entered as splines, and β q (q, x) or β r 1 r 2 (r 1 , r 2 , x) functions modeled as two-or three-dimensional tensor products, respectively. Spatial coordinates were projected and included in the model to account for spatial autocorrelation and to model geographic patterns not represented in other covariates. The maximum number of knots that could be used for estimating spline functional terms was set at 11 as a balance between overfitting and modeling abrupt changes in non-crop area over annulus distance. All covariates were scaled and centered.
Interpretation of the model involved examining the coefficient functions α p(i) (x) and β q q, x or generating predictions by excluding certain terms and systematically selecting input values to examine scenarios of interest. For example, predictions controlling for environmental covariates were generated by setting the relevant categorical functional or continuous functional terms to the zero function and entering categorical parametric coefficients as constants at the level with the most observations. Mapping was done by including the spatial term, β r 1 r 2 (r 1 , r 2 , x), but not other environmental covariates, and generating predictions for the centroids of survey townships. The potential for increasing landscape complexity was assessed by finding the difference in the mean non-crop area prediction between neighborhood and within-field scales for the centroid of each survey township.

RESULTS
The data set represented 14,527 randomly-selected quarter sections with a mean crop area of 79%, and a mean distance from nearest neighbor of 2,891 m. These quarter sections occupied 2,430 survey townships with a mean of 6 quarter sections in each (min = 1, max = 14), collectively sampling from an agricultural study region of 226,543 km 2 . Forty-five annuli ranging from 150 m to 1,006 m in radius were assessed for each quarter (mean area per annulus = 7.1 ha). Data from the sampled quarters and their surroundings covered 6% of the total area under study.
The model deviance explained was 32.2% (R 2 adj = 0.355) with all constant parametric coefficients for categorical terms significantly different from zero (P < 0.01; Table 2a) suggesting the average proportion of non-crop area differed by ecoregion and agricultural limitation class. All functional non-parametric coefficient functions significantly differed from zero (P < 0.01; Table 2b) demonstrating a non-linear relationship between noncrop area and the distance from field centers, and that the included covariates mediate this relationship. The spatial term summarizing the effect of unmodelled geographic variables also had an effect on the landscape complexity curve (Northing × Easting; Table 2b).
The average trend in the proximity of non-crop area to field centers is that inside fields (150-400 m) there tend to be lower proportions of non-crop area than in the surrounding landscape (500-1,006 m; Figure 3). Thus, focal crop fields were not typically surrounded on all sides by other crop fields (a situation that would produce no difference between these two scales). However, there is spatial variation across the region in the proximity of non-crop areas to fields, after controlling for environmental covariates. This pattern is evident in the geographic differences in the predicted amount of non-crop cover at the withinfield scale (e.g., Figure 4A), and at the neighborhood scale (e.g., Figures 4C,D). The field margin scale (400-500 m; e.g., Figure 4B) is more conserved across the province.
Most ecoregions ( Figure 5) and agricultural limitation classes (Figure 6) had a characteristic landscape complexity curve, with FIGURE 3 | The model intercept function (or mean landscape complexity curve) shown with two standard errors. The curve represents the estimated effect on the proportion of non-crop area and how it varies over distance from the center of an average field. Vertical dotted lines in this and subsequent figures indicate annulus radii at which fields were measured, providing support for the model.  these coefficient function plots demonstrating deviations in the coefficient from the mean landscape complexity curve (i.e., from Figure 1). All climatic variables also exerted effects on the shape of the landscape complexity curves (Figure 7). Trends for scales larger than the field margin scale (400-1,006 m) were generally similar for temperature and moisture-related variables. At the within-field scale, higher temperatures, and drier soils resulted in less non-crop area (mean annual temperature, plant available water; Figure 7).
Mapping of a multi-scale index shows that there is geographic variation in the difference between the proportion of non-crop covers at the neighborhood and the within-field scales. Larger

DISCUSSION
This study demonstrated that there is considerable geographic variation in the proximity of non-crop areas to fields across Alberta's agricultural zone, and there remains capacity for growers and other land use decision makers to optimize this distribution. Analysis of landscape complexity curves revealed scale-dependent patterns in the proximity of non-crop areas to fields, such as the widespread occurrence of uncultivated field margins, and how land uses have been influenced by broad environmental gradients. These findings are considered in turn.

Potential for Changing Landscape Complexity
The capacity to introduce more non-crop land covers into fields, and therefore improve landscape complexity, can be seen by contrasting the proportions of non-crop area found within-fields to the neighborhood surrounding those fields (Figure 8). Noncrop area at the neighborhood scale can be understood as an estimate of the local potential for this quantity. Annuli at this scale sample from eight neighboring quarter sections, many of which may also be fields, effectively summarizing a region eight times the size of the focal field (Figure 1). The neighborhood scale, then, can be taken as a realistic potential proportion of non-crop area in the focal field because it captures neighboring land owners' willingness to allow those non-crop land covers to persist. It is also likely to be a correlate of the local expectation for crop productivity, given that neighboring growers would have an incentive to drain wetlands, clear trees and shrubs, remove perennial grasses, cultivate to fence lines or otherwise remove non-crop land covers from their fields, if that land could be used profitably for growing crops.
For example, the model suggests that within fields at location P (annulus = 150 m; Figures 4A,E) there is about 20% coverage in non-crop area, and there is also a similar coverage in the neighborhood (annulus = 995 m; Figures 4D,E). While 20% coverage suggests that non-crop covers are relatively common in fields, the more interesting observation is that there tend to be adjacent fields with similar non-crop proportions. This emerges as a flatter curve, and it is interpreted here as evidence that landscape complexity has reached its local potential; i.e., that land owners generally concur that 80% of the land should be in crops. The contrasting trend is evident at location T, where there is lower coverage of non-crop area within focal fields (about 5%) and higher coverage at the neighborhood scale (about 25%). Here fields have a low amount of non-crop area, but there is a large difference between the field and the neighborhood. This appears as a steeper curve, offering evidence that fields have less non-crop area than the local potential; i.e., that land owners differ markedly in how much land should be in crops. Figure 8 (darker shades) therefore identifies parts of Alberta where agricultural land uses may be out of equilibrium with what the landscape can sustain, either because there may be more of each field in crop production than is optimal, or because there remains potential to further clear natural land covers for agriculture. The latter is probably the case in northern areas where twenty-first century clearing of boreal forest for FIGURE 8 | The difference in the predicted non-crop area between the within-field (150-400 m) and neighborhood (500-1,006 m) scales for each of the 2,480 survey townships with agricultural land. Model predictions exclude all environmental covariates but include spatial terms.
agriculture is ongoing (north of 55 • N; Bowen, 2002). However, the opportunity to clear more land is unlikely to be broadly important further south where agriculture has had an impact on the land since the late nineteenth century (Larmour, 2005).
The parts of the province that have the greatest difference at these two scales are those between 52 • N and 55 • N, after excluding the northern region and a region in south-central Alberta that is predominantly forage land (e.g., near point T; Figure 4). These latitudes are mostly in the parkland ecoregions, where prime agricultural land is abundant (Figure 2). There, growers may be able to leverage the collective wisdom of their neighbors to guide the return of non-crop areas to their fields. Or, in more concrete terms: marginal areas of fields that are currently cropped could be changed to non-crop vegetation; and forgoing production in these areas would mirror the decisions that neighbors have, on average, made regarding land use. This study, however, cannot provide any insight into how much noncrop land cover could be recovered. Rather it demonstrates that there are parts of the province where patterns in land use suggest the barriers doing so should be smaller than others.
Making the assumption that non-crop land covers in fields improve ecosystem services to agriculture (e.g., Rusch et al., 2016;Duarte et al., 2018;Vickruck et al., 2019), there is also the possibility that a small loss in crop area and therefore in farm returns, may be balanced by the improvement in the yields and profitability of crops on the remaining land. This can further reduce the risk associated with removing crop and replacing it with non-crop land covers.
Increasing non-crop land covers within fields could, in many parts of this region, be achieved by identifying marginal areas with relatively low productivity (e.g., measured in crop yield), removing them from production, and allowing other noncrop vegetation to re-establish. The proliferation of precision agricultural tools (e.g., yield monitors in harvesting equipment) should assist growers in identifying such sites (Mulla, 2013). These may be low spots in fields that have excess soil moisture or other poor soil conditions, or they may be near to other features that reduce productivity such as in the margins of wetlands where soils are poorly optimized for crop growth.
This study has not given any consideration to the class of noncrop cover and its association with landscape complexity, in part because this simplification aids interpretation. However, types of vegetation may differ in the ecosystem services they can provide to nearby crops, suggesting that regionally-targeted research on the benefits of establishing different vegetation classes is certainly appropriate. Equally, biasing re-establishment efforts toward perennial plants may better support carbon storage objectives (Hungate et al., 2017;Williams et al., 2018). However, the classes of vegetation already common at the neighborhood scale may be those that are the easiest to re-establish, and these may also be the species that naturally reclaim these sites without any intervention from land use decision-makers.

Patterns in the Proximity of Non-crop Area to Crops
Analysis of landscape complexity curves reveal field margins (400-500 m from field centers) as a common location for noncrop cover. These are evident as a peak at locations P, Q, R, S, and T ( Figure 4E) at the scale that corresponds to the expected survey grid spacing. Geographically, the field margin effect is found throughout much of the region, with the notable exception of the extreme south of Alberta (annulus = 474 m; Figure 4B).
Field margins have been celebrated as a means to maintain non-crop areas in agricultural landscapes and to bring the ecosystem services they may provide close to fields with minimum loss of crop area (e.g., Marshall, 2002). In many cases these field margins are adjacent to roads or road allowances, which are public lands and are therefore at low risk of being changed to other land uses. Their widespread geographical distribution in Alberta represents an opportunity for regional policies that systematically promote their enhancement (e.g., by altering mowing regimes, or planting with native vegetation known to support ecosystem service provision; Kirmer et al., 2018).

Environmental Drivers of Landscape Complexity
Broad environmental gradients have played a role in shaping landscape complexity resulting in different patterns of non-crop area across scales. As might be expected, there is evidence that these gradients have influenced land use decisions, particularly the amount of clearing within fields and the density of fields at the neighborhood scale. Boreal ecoregions tended to have more noncrop area within fields (150-400 m; Figure 5) than the regional mean (Figure 3), perhaps reflecting the more recent clearing of these northern areas (Bowen, 2002) and the greater speed with which shrubs and trees can re-establish in the moisture and temperature regime of this part of the province. Grassland and parkland ecoregions registered at or less than the mean at the within field scale. The Mixed Grassland ecoregion had the most vegetation in field margins (400-500 m) while two ecoregions with significant tree cover (Mid-Boreal Uplands, and Boreal Transition) had the most non-crop area at neighborhood scales, reflecting the substantially lower amount of agricultural activity in these regions (Figure 5; Ecological Stratification Working Group, 1995).
Agricultural limitation classes demonstrated there has been more removal of vegetation within fields where there are fewer challenges to crop production in terms of the workability of the soil (Figure 6). Crop fields that have been successfully established in areas with severe limitations to agriculture suggest they are similar to the surrounding landscape in terms of non-crop area, and field margins are not distinguishable (e.g., marginal agriculture; Figure 6).
Temperature, precipitation and soil moisture variables enabled a directional assessment of how climatic variation is associated with landscape complexity patterns. Notably, locations with warmer temperatures and less soil moisture, conditions which favor prairie grassland species (Ecological Stratification Working Group, 1995), had less non-crop area within fields (Figure 7). The low frequency of woody vegetation in such dry prairie conditions means that clearing fields of all non-crop area is easier and less costly. Overall, the studies of environmental gradients indicate that landscape complexity is primarily under the control of land use decision-makers and not merely a result of local conditions, suggesting it is a matter of incentivizing these changes and not a problem simply of what the environment can sustain.
The multi-scale approach used in this study provided a flexible means to characterize landscape complexity as the proximity of non-crop features to crop fields. Interpretation of this model revealed that there is potential to increase landscape complexity by leveraging the natural potential in each region to support non-crop land cover and promoting practices that foster such vegetation, for example, on marginal or low-productivity sites within fields.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
PG and MG were involved in the writing of the manuscript. PG was primarily responsible for the completion and interpretation of analyses. MG compiled the composite land cover data product on which this analysis was based.