**Leveraging Traditional Agroforestry Practices to Support Sustainable and Agrobiodiverse Landscapes in Southern Brazil**

#### **André Eduardo Biscaia Lacerda 1, Ana Lúcia Hanisch <sup>2</sup> and Evelyn Roberta Nimmo 3,\***


Received: 11 March 2020; Accepted: 28 May 2020; Published: 1 June 2020

**Abstract:** Integrated landscape approaches have been identified as key to addressing competing social, ecological, economic, and political contexts and needs in landscapes as a means to improve and preserve agrobiodiversity. Despite the consistent calls to integrate traditional and local knowledge and a range of stakeholders in the process of developing integrated landscape approaches, there continues to be a disconnect between international agreements, national policies, and local grassroots initiatives. This case study explores an approach to address such challenges through true transdisciplinary and multi-stakeholder research and outreach to develop solutions for integrated landscapes that value and include the experience and knowledge of local communities and farmers. Working collaboratively with small-scale agroforestry farmers in Southern Brazil who continue to use traditional agroecological practices to produce erva-mate (*Ilex paraguariensis*), our transdisciplinary team is working to collect oral histories, document local ecological knowledge, and support farmer-led initiatives to address a range of issues, including profitability, productivity, and legal restrictions on forest use. By leveraging the knowledge across our network, we are developing and testing models to optimize and scale-out agroforestry and silvopastoral systems based on our partners' traditional practices, while also supporting the implementation of approaches that expand forest cover, increase biodiversity, protect and improve ecosystem services, and diversify the agricultural landscape. In so doing, we are developing a strong evidence base that can begin to challenge current environmental policies and commonly held misconceptions that threaten the continuation of traditional agroforestry practices, while also offering locally adapted and realistic models that can be used to diversify the agricultural landscape in Southern Brazil.

**Keywords:** yerba mate; agroforestry; integrated landscape; agrobiodiversity; silvopastoral systems

#### **1. Introduction**

Several high-level reports from a range of international agencies highlight the need to rethink conventional agricultural systems through innovative and sustainable approaches, including agroecology, forest landscape restoration, and agroforestry, among many others [1–6]. These reports emphasize that business as usual in terms of conventional agriculture is continuing to have lasting negative impacts on agricultural biodiversity, soil health, water and landscape management, greenhouse gas (GHG) emissions, and human health and food security. The United Nation Food and Agricultural Organization's (FAO) recent report on the state of the world's biodiversity for food and agriculture (BFA) argues that "many of the drivers that have negative impacts on BFA, including overexploitation, overharvesting, pollution, overuse of external inputs, and changes in land and water management, are at least partially caused by inappropriate agricultural practices" [2] (p. xxxviii). Meanwhile, Padoch and Sunderland [7] argue that using conventional practices and technologies for sustainable intensification may not necessarily have the desired effects on forest and biodiversity conservation, but rather lead to greater loss of forests and associated ecosystem services, with little or no benefits for some agricultural regions in which small-scale farming is predominant. They highlight that "producing food in diverse, multifunctional landscapes challenges dominant agricultural development paradigms, but it also presents issues and difficulties. For example, many types of integrated landscape approach have not been studied by scientists, and the existing research and policy framework may be insufficiently integrated to improve either agricultural production or environmental protection in such diverse landscapes" [7] (p. 6).

From a landscape ecology perspective, taking a holistic approach to land management planning and modeling is a key aspect of the discipline [8,9]. Understanding the mechanisms and impacts of land use and land cover change (LULCC) over time in farming regions, including forest fragmentation, habitat loss, and human–environment interactions, is crucial to determining the likely impacts of the continuation of conventional agriculture not only at the local scale, but also how this will affect rural and urban human and non-human populations in the long term. In order to support sustainable management and changes to land use and land cover (LULC) that focus on biodiversity conservation, increased ecosystem services and connectivity, as well as human food security and livelihoods, debates have focused on land sharing vs. land sparing as methods to address the competing needs in landscapes [7]. This debate either calls for land to be set aside for conservation with intensive agriculture conducted separately, or land to be shared across a range of goals from food production to biodiversity conservation through less intensive practices such as agroecology [10]. Although many involved in the debate acknowledge that both land set aside for conservation and alternative approaches to agriculture can occur simultaneously, there are few examples of how this might work on a practical level or how to scale up what works on individual farms to address issues of managing sustainable, biodiverse productive landscapes.

One of the methods used in landscape management planning that bridges the divide across the various competing social, ecological, economic, and political contexts and needs in landscapes are integrated landscape approaches. As defined by the Consultative Group on International Agricultural Research (CGIAR) [11], such approaches consider not only multiple land uses (including agriculture and forests) but also the livelihoods dependent on such land uses, moving beyond conventional perceptions of management and governance. It seeks "to provide tools and concepts to identify, understand and address a complex set of environmental, social and political challenges, and to enable evidence-based and inclusive prioritization, decision-making and implementation" [6] (p. 1). Importantly, such an approach highlights stakeholder engagement as key to managing conflicting perceptions of the value and function of land use types in a landscape across a range of scales, from the local to the national [6]. What is important here is that analyzing and developing solutions for integrated landscapes requires a truly transdisciplinary lens, in which a range of researchers and other stakeholders, including local communities and farmers, are actively engaged in the research design, data collection and analysis, implementation, and assessment.

Agroforestry systems have been identified as one of the key approaches that can be implemented in integrated landscape management as they offer a range of ecological and social benefits. As noted in the recent High Level Panel of Experts (HLPE) report, forests contribute extensively to food and nutritional security (FSN) through the "direct provision of food; provision of energy, especially for cooking; income generation and employment; and provision of ecosystem services that are essential for FSN, human health and well-being" [6]. Specifically, the implementation of agroforestry systems offers multifunctional landscapes that support the development of regenerative agricultural systems that offer not only a diverse, multi-layer food production system, but also land use that can restore or conserve ecological resources [12]. The International Union for Conservation of Nature (IUCN) together with the World Resources Institute (WRI) have highlighted agroforestry methods as an

important strategy in forest landscape restoration (FLR) to address climate change mitigation and food security issues worldwide [4]. Research on the benefits of agroforestry systems and their ecological, social, economic, and cultural importance have been conducted in a range of contexts around the world (see for example [13,14]).

In Southern Brazil, previous research has shown that forest fragments, including those managed in agroforestry systems, are important havens of biodiversity on the landscape scale, particularly in terms of tree species [15–18], but also act as crucial connectivity corridors that enable genetic flows and buffer the impacts of anthropogenic activities along waterways [19]. Traditional agroforestry systems have continued in the central-south of Paraná state and Northern Santa Catarina state mainly due to the extraction of erva-mate (*Ilex paraguariensis*, also known as yerba mate), a tree species that grows well in the shaded understory of the region's iconic Araucaria Forest. These systems have been important in maintaining ecosystem services and biodiversity corridors, but they are also important to the maintenance of cultural and traditional agroecological practices on small-scale family farms that include a heterogenic mosaic of crops, livestock, vegetable gardens, and productive forest areas, all of which are essential to family and local food security [20].

We understand these traditional systems as outcomes of generations of adaptive practices based on local ecological knowledge (LEK), resource management techniques, and cultural and historic subjectivities [21]. We leverage the premise outlined by Berkes et al. [22] and Fonseca-Cepeda et al. [21] of traditional ecological knowledge (TEK), which Berkes [23] (p. 3) defined as "a cumulative body of knowledge and beliefs, handed down through generations by cultural transmission, about the relationship of living beings (including humans) with one another and with their environment. Further, TEK is an attribute of societies with historical continuity in resource use practices." Although such a concept is often associated with indigenous knowledge paradigms, we consider how such an approach can apply to settler communities that have continual, historical resource use practices, as this knowledge is "cumulative and dynamic, builds on the experience of generations", but also adapts to new technological and socioeconomic realities [24] (p. 281). Nevertheless, these systems have received little research attention, particularly from federal and state agricultural research and outreach institutions, as they tend to be viewed as remnants of outdated, subsistence agricultural practices that require modernization, rather than being valued as systems developed and adapted to the forest environment in which they have continued for generations.

In this paper, we discuss some preliminary results of our ongoing participatory research and outreach project with small-scale traditional agroforestry producers in Southern Brazil and present models being developed to optimize and scale-out agroforestry systems based on our partners' traditional practices. Through the implementation of these models and systems, we are beginning to address some of the main concerns of small-scale farmers, mainly profitability, productivity, and legal restrictions on forest use, while also developing strategies that can be used in landscape management to expand forest cover, increase biodiversity, protect and improve ecosystem services, and diversify the agricultural landscape. Our collaborative approach ensures that the research addresses the needs of communities and is applicable to local realities. In so doing, we are developing a strong evidence base that can begin to challenge current environmental policies and commonly held misconceptions that threaten the continuation of traditional agroforestry practices.

#### **2. Traditional Land Use and Conventional Agriculture**

The land use legacy of Paraná state helps to clarify the current LULC in the region and how the continuation of forest resources in some regions has been more pronounced than others, as shown in Figure 1. At the beginning of the nineteenth century, most of Southern Brazil was covered by forests, from the coastal Atlantic Forest through to the sub-tropical Araucaria Forest biome on the central plateaus, and the tropical semi-deciduous forests of the Paraná River basin in the west. Although the forests had been managed by indigenous groups for thousands of years [25] and while there was continued indigenous and settler occupation of these forest landscapes from the sixteenth century [26], forest cover was relatively uniform, interspersed with the naturally occurring grasslands on the higher elevation plains. By the end of the eighteenth century, a process of westward colonization began in Paraná, and to some extent Santa Catarina state, originating from the coastal/eastern region and moving through the region's highlands. This colonization process was characterized by an economy based on cattle husbandry, erva-mate harvesting (a resource that had only recently begun to be exploited in the regional economy, despite its economic importance since early Spanish colonization in the seventeenth century [27]) and the logging of araucaria or Paraná pine (*Araucaria angustifolia*) [28]. More intense colonization did not happen until later, with a migratory influx of Germans, Italians, Poles, and Ukrainians, among others, in the nineteenth and early twentieth centuries, as part of a government policy to occupy the 'unoccupied' hinterlands of the country [29]. This second wave of migration, along with the expansion of the railroad and extensive exploitation of araucaria, led to much more intensive land occupation. Yet, in our study region, many settler communities maintained forest cover on their lands for animal husbandry and erva-mate harvesting, developing farming systems that integrated European and indigenous practices with local crops (corn, manioc, beans, etc.) and forests. Today, forests still occupy significant portions of the landscape, as shown in Figure 1c, in which small-scale farmers continue to use traditional practices that have been passed on for generations.

**Figure 1.** Forest and land use and land cover (LULC) in Paraná state as a result of differing land use legacies: (**a**) location of the region in Brazil; (**b**) Paraná state with northern and southern regions highlighted; (**c**) land use in southern region with a significant incidence of forest cover; (**d**) land use in northern region mostly covered by soy plantations.

In other areas of Paraná, where forest cover has been almost completely decimated, the process of colonization was much later and with quite a different focus. Beginning in the early to mid-twentieth century, the Southwest was occupied by colonization companies expanding from the southernmost state of Rio Grande do Sul; their economies were linked to a subsistence agriculture based on grains and pork husbandry. Finally, a colonization wave coming from the north (São Paulo state) in search of land for coffee production occupied the north of Paraná [29]. The vast majority of the area originally occupied by the colonization from São Paulo and Rio Grande do Sul is currently covered by large-scale soy farms, as shown in Figure 1d.

#### **3. Methods to Leverage Traditional Agroforestry Practices**

Our approach to implementing participatory methods to develop land management planning systems that integrate multiple uses of natural resources in a socially and politically complex context was conceptualized as Locally Adapted Participatory Sustainable Forest Management Systems—lapSFM, outlined by Lacerda et. Al [30]. This approach provides a roadmap for managing rural properties, focusing on forested lands, while the decision-making process includes local ecological knowledge (LEK) as a key input necessary for establishing the goals and objectives and is based on the demands and interests of landowners. As such, in 2011, we began by leveraging long-term agroecological initiatives in place for more than 30 years that involved key partners (Federation of Family Farmers' Unions—FETRAF; Agronomic Institute of Paraná—IAPAR; and the Brazilian Agricultural Research Company—EMBRAPA), and started to conceptualize these productive systems by consolidating communication, information, and activities among farmers, outreach technicians, researchers, and environmental agencies through workshops, field-days, social media, farmers' union meetings, and scientific conferences.

Beginning in 2017, we began a new phase of the project which focused on conducting oral history interviews with traditional erva-mate producers as a means of gaining a better understanding of the historical, social, and cultural aspects of traditional practices, how landscape and environmental changes are perceived by farmers, as well as the socio-political implications of their lived experiences [31]. As Williams and Riley [32] argue, oral history interviews provide an understanding of the ways in which people produce meaning of the places they inhabit, and how they perceive and value the natural world around them. As such, they offer unique perspectives on issues of the environment, forests, and conservation as narratives are situated within the environment in question, grounded in the everyday challenges of rural life. To date, we have conducted interviews with 39 erva-mate producers and members of their families across seven different municipalities in Paraná and Santa Catarina [33].

Across the range of participatory methods employed, participants were encouraged to discuss the challenges faced in conducting agroecological and traditional agroforestry practices on small-scale farms (economic, technical, social, and political) and possible solutions that included creation of co-ops, cooperative/participatory research to deal with gaps in scientific knowledge, youth-focused farming, and increased participation of women, among others. As a result, we have established action plans and models that integrate a range of perspectives and issues in terms of the technical, social, cultural, and economic.

Herein, we discuss two particular outcomes of project to date: optimized LEK-based agroforestry systems; and models of Productive Agroforestry Restoration. These optimized productive systems have been implemented at EMBRAPA Research Station in Caçador (ERSC) and in over 50 properties across more than 20 municipalities in Paraná, Santa Catarina, and Rio Grande do Sul states, for a total of 3000 hectares under management. Both lines of research seek to add value to these systems and the knowledge behind them, while also testing alternatives that do not require high rates of investment or debt, and may offer small-scale farmers the opportunity to transition some of their land from high-input commodities (i.e., corn, tobacco, soy) to other more sustainable products. The dissemination of results and co-creation of knowledge include technical and scientific documents [19,20,34–40], conferences e.g., [41], monthly technical visits to ERSC, and bi-monthly visits to farms across the region.

#### *3.1. Traditional Agroforestry Optimization*

The diversity of forest management systems based on the traditional use of erva-mate reflect the variations in the natural environmental that led to an extensive accumulation of LEK over generations. Forest structure, tree diversity, presence of dominant or invasive species, and land use history and legacy are all integral factors that play a role in the decision-making process related to how forests are managed for erva-mate production and have been described in depth by Mattos [42], Chaimsohn and Souza [18], Marques [43], and Hanisch [44]. In most cases, forest structure and diversity are gradually managed, aiming at a spatial distribution that favors an understory with a homogeneous light availability considered empirically as ideal for erva-mate development. The intensity and frequency of forest interventions depend on various factors that include forest development (successional stage), historic and contemporary use, presence of dominant or invasive species, and natural occurrence of erva-mate trees, among others. Despite the fact that traditional agroforestry systems in Southern Brazil in most cases have erva-mate as one of the key products, practices vary widely between properties and municipalities, with the presence of cattle husbandry as one of the most significant characteristics that differentiates the systems in the region [20].

#### 3.1.1. Agrisilvicultural Systems

Erva-mate production occurs across a range of forest successional stages, from well-developed (late successional) forest stands relying mainly on native, naturally regenerating trees, to younger, secondary forests that rely more heavily on planting and silvicultural management. Well-developed forest stands typically have lower density (number of trees), higher diversity with long-living species (i.e., *Araucaria angustifolia* and *Ocotea porosa*), with a more complex structure (trees with various sizes distributed in forest layers), while younger forests commonly have a much simpler structure (homogeneous sizes), lower levels of diversity with short-lived species (i.e., *Mimosa scabrella*) but with much greater density. Research has shown that agroforestry systems with erva-mate in Southern Brazil present significant levels of tree diversity, with 107 tree species identified across 39 botanical families [18], which represent a significant proportion of the region's diversity [15]. They also show high levels of nutrient cycling through litter that far exceed the nutrient exportation that takes place during erva-mate harvesting [45].

In well-developed forests, erva-mate occurs naturally as large trees and is harvested through radical pruning in 2- to 4-year production cycles. Erva-mate can be also be planted and managed at shrub size for ease of harvesting, but production is often limited due to low levels of luminosity for plant development. Although well-developed forests are mostly found in an "open" state [43], in which historic management has reduced tree density, as shown in Figure 2a, some management is required to ensure optimal conditions for erva-mate growth. However, canopy management as a means to increase light in the understory is limited by very restrictive legislation governing the use of such forests (see the Brazilian Forest Code [46,47] and Atlantic Forest Law [48]). Consequently, production in well-developed forests tends to be restricted to animal husbandry with low stock density and naturally growing erva-mate plants. Contrastingly, younger forests tend to be more actively managed in Southern Brazil for erva-mate production. In most cases, the understory is maintained mostly clear in order to make space for erva-mate plants, while thinning is applied if insufficient light permeates the forest canopy. The intensity of thinning is established empirically and ranges widely from intense tree reduction where producers aim at production levels similar to monoculture stands, to agroecological practices that try to maintain sustainable multispecies environments.

**Figure 2.** Traditional agroforestry optimization: (**a**) traditional erva-mate production where the lack of adequate forest management caused canopy decline leading to monoculture-like erva-mate production; (**b**) understory management, with initial removal of invasive native bamboos, for production and species diversification; (**c**) induced forest regeneration for canopy restoration (white trunks), followed by erva-mate plantation (increased production) and diversified tree plantation (species and forest structure diversification); (**d**) sheep used for weed control in a newly implemented erva-mate intensification stand.

Our research on the management of young forests has led to the establishment of best practices that aim to achieve a more stable income while maintaining or increasing ecosystem diversity and complexity using innovative silvicultural treatments focused both on erva-mate plants and management of other forest elements [39]. In this method, erva-mate is harvested annually at much lower intensities (~50%) that maintain plant vigor as opposed to the 2- to 4-year cycles during which almost all plant foliage is removed, causing significant plant stress. Site-specific silviculture is used to define the plant's shape, height, and pruning, and regeneration tending. We introduced the use of spreaders (strings, bamboo) to widen erva-mate tree crowns into a goblet shape that allows for a greater production per plant compared to traditional growth based on a cluster of vertical branches growing from a main trunk (lower leaf-to-branch ratio). Furthermore, branches are trained away from the main vertical axis of the tree creating a sharp to near horizontal angle from which new sprouts can be harvested. Leaf harvesting from secondary branches instead of a radical pruning from the main trunk dramatically reduces the damage caused by water seepage that rots the trunk interior and often occurs after consecutive pruning, ultimately reducing the plant's life-span.

Additionally, areas where light is more available at the understory level will have plants pruned and harvested at 1–3m in height, whereas a darker understory will require taller trees (3–8m) to improve production. Similarly, spacing between trees will determine pruning methods: sparser plantations (e.g., 3 × 2 m) allows for shaping a much wider crown which in turn produces abundant foliage that can be harvested annually at less stressful levels (~50%).

We also introduced practices to ensure forest regeneration for maintaining forest cover and diversity in the long term, avoiding the expensive and impractical need to reintroduce trees through seedling planting. Two simple methods were successfully tested [39] and currently applied in farms: (i) identification and marking of seedlings from natural regeneration using bamboo/wood sticks prior to weed trimming; (ii) defining areas in which weed trimming is not conducted—plots of 1 m2 are marked using sticks where natural regeneration is protected and encouraged. Our results show that after one year, at least one native tree seedling was successfully recruited in 75% of these plots (reaching up to more than 20 recruits in a single plot). We recommend that such areas of regeneration are introduced together with the planting of erva-mate, replacing one erva-mate seedling every 5–9 m to support forest succession.

One of the opportunities identified through this research and by others [43] is the management of young forests in the region that are dominated by native invasive bamboo species, as shown in Figure 2b. Bamboos are a determinant factor in forest dynamics as they tend to impede the development of seedlings and young trees [37,49,50]. Our previous research has shown that these bamboo species, when dominant, create a simplified plant community in which succession is arrested [37,51] and an impractical environment for the development of productive systems. In terms of landscape management planning, the areas in which bamboo are dominant require practical management solutions that kick-start forest succession to improve biodiversity and maintain forest connectivity, but also offer benefits to property owners. In these areas, farmers can manage the understory, eliminating bamboo cover and establishing erva-mate plantations in densities varying from 1500 to 10,000 seedlings per hectare, as shown in Figure 2c. Again, farmer's objectives define specific practices: production maximization tends to lead to very dense plantations with the use of chemical fertilizers and pesticides (an illegal process as no pesticide is regulated for use with erva-mate) along with the removal of forest regeneration and continuous canopy thinning; whereas traditional and agroecological producers use organic (or no) fertilizers with forest management aiming at maintaining forest structure and diversity in the long term. Typically, the former produces more leaves per area, while the latter often obtains a premium for the quality of the product.

New techniques to improve traditional agroforestry practices are also incorporating farmer-led initiatives that have found innovative ways to optimize production and minimize environmental impacts. One very promising solution is the use of sheep husbandry for weed control in erva-mate plantations, as shown in Figure 2d. The introduction of the sheep breed "Texel" for controlling weeds reduces substantially the demand for labor, which is one of the most pressing limitations in farming today. One important characteristic of the Texel breed is the fact that they do not graze on erva-mate plants and have minimal impact on soil compaction. Finally, sheep farming can play an important role in food security for farmers and can be a smart solution to halt the current trend of using glyphosate herbicides for weed control. Such low-tech, practical approaches to optimizing production and reducing the use of chemical inputs can be scaled out from individual small-scale farms to create regional approaches and best practices that recognize and support the knowledge and participation of farmers, and in turn can have lasting impact on the forested landscape.

#### 3.1.2. Silvopastoral Systems

The use of animal husbandry in traditional agroforestry in Southern Brazil has two main systems: *caívas* in North Santa Catarina state and *faxinais* in South Paraná state. *Faxinais* were once a common feature in the landscape in Paraná, where communities use a large forest area as a commons for animal husbandry, as shown in Figure 3a, and erva-mate harvesting, with a wide diversity of food crops (including corn, beans, manioc, and rice) in fenced fields protected from animal grazing [43]. These multifunctional land use systems were typical of the indigenous descendants that occupied the region and later assimilated by settler communities, especially Ukrainian and Polish immigrants. On the

other hand, *caívas* are generally found in the Northern Plateau region of Santa Catarina on individual properties in which dairy cattle husbandry is carried out in forest patches usually combined with erva-mate production, as shown in Figure 3b [20].

**Figure 3.** Traditional silvopastoral systems: (**a**) animal husbandry within forest stands in a *faxinal*; (**b**) *caíva* with low productivity pasture and senile erva-mate trees (small trees in the background, to the right); (**c**) view of a *caíva* with improved pasture.

The natural pasture in *caívas* typically has low levels of productivity, particularly in the winter when plant regrowth cannot keep up with grazing demands [44,52]. Thus, animal productivity is low, leading to food insecurity and low economic income and resilience. Undernourished animals have knock-on economic impacts on erva-mate production as they look for grazing alternatives and damage erva-mate trees or consume the leaves. Furthermore, environmental sustainability may also be affected because grazing on forest regeneration can compromise forest renovation and physical damage caused by bark consumption, which compromises tree health, ultimately reducing tree lifespan.

As farmers look for more profitable economic alternatives, *caívas* have been replaced by monoculture production based on forest plantations and commodity crops with direct loss of forest cover and biodiversity and traditional practices. As a response, a participatory research project carried out by EPAGRI, the Agricultural Research and Rural Outreach Company of Santa Catarina, and EMBRAPA has focused on developing strategies to increase household income by improving animal and erva-mate productivity while also maintaining or restoring forest structure and diversity [20,44,52]. The project framework includes testing innovative practices related to pasture improvement, renovation of forest stands and erva-mate trees, and defining ideal levels of canopy cover and regeneration management.

In 2010, EPAGRI initiated the implementation of improved practices for traditional *caívas* that include pasture overseeding during the winter that evolved into the development of the genetically improved *Axonopus catharinensis* SCS 315 (referred to locally as *missioneira-gigante* or giant missionary grass) [53], a pasture species that is better suited for *caíva* environments due to its tolerance to shade [44,54,55]. The process to implement this technique includes the removal of native grasses and the introduction of giant missionary grass (detailed information about the process can be found in [44]; Figure 3c). The improved *caíva* system has been implemented successfully in eight farms across four municipalities in Santa Catarina [44] and another ten properties will adopt the technology by the end of 2020; those farms will be used as reference properties for outreach agencies to disseminate the results.

Along with the goal of increasing animal productivity, the project also developed practices to improve erva-mate production. These include a set of activities aimed at creating a highly diverse, healthy forest with multi-aged and multi-strata elements. Due to years of neglect and prohibitive laws that severely restrict forest management, as noted above, most traditional *caívas* have inconsistent forest structure ranging from large, frequent gaps to very dense clusters of trees. Thus, we developed forest management guidelines to support forest restoration, canopy refinement, and erva-mate intensification. Forest restoration seeks to restore gaps in forest cover and reintroduce a multi-aged tree population which can be achieved by designating areas for restoration where animal grazing is temporarily restricted (usually by using electrical fencing) for 3–5 years, after which foraging is again allowed and a new area is fenced. Monitoring of these areas showed that regeneration was highly effective in restoring species diversity and structure as 59 different tree species were recorded in the fenced areas [56,57], which was greater than the diversity of the adult tree population. Additionally, we recommend thinning of abundant species and tree clusters in order to increase species diversity and establish a more even forest canopy, respectively. Simultaneously, in fenced areas, erva-mate can be planted at densities between 1000 and 3000 seedlings per hectare to increase productivity.

Through the implementation of the higher-productivity, shade-adapted perennial grass in *caívas*, farmers are producing five times more pasture per area, enabling a triplication in the stocking rate and consequent increases in milk production [43]. Furthermore, erva-mate production can increase tenfold depending of the level of degradation of the erva-mate trees. Finally, in order to monitor changes over time and help evaluate the impacts of new practices as they are implemented, we adapted the Sustainability Assessment of Food and Agriculture Systems (SAFA) developed by FAO [58] to assess the sustainability of farms considering 77 indicators across themes of environmental integrity, economic resilience, quality of life (social), and good governance [20]. The overall results showed that the implementation of the improved practices outlined above enabled farmers to obtain a ranking of good (based on a classification as unacceptable, limited, or good) for 87% of indicators, in comparison to 65% ranked as good for *caívas* that had not implemented improved practices.

As is the case in many other countries [7], current agricultural and environmental policies in Southern Brazil do not recognize traditional practices, as they are often excluded from scientific analyses or assessments and intensification through conventional agriculture is still seen as the way forward. As such, landscape policies and management strategies do not consider how traditional systems might be leveraged to mitigate the increasing homogeneity of the landscape through monoculture farming and the consequential impacts on human health and nutrition, biodiversity, gene flow, and ecosystem services. Without a recognition of, or support for, traditional systems, the tendency is for this local knowledge to be lost, along with the associated cultural identities and environmental subjectivities. Thus, analyses, participatory approaches, and farmer-led initiatives such as those we have outlined in this section that optimize traditional systems, provide the evidence base necessary for these practices to be integrated into policy and governance structures, which in turn can provide landscape managers with practical approaches that are culturally relevant and can be implemented across the landscape.

#### *3.2. Productive Agroforestry Restoration*

As noted above, Southern Brazil has been subjected to devastating rates of deforestation and forest fragmentation over the last century, as shown in Figure 1. While current policies have been essential in reducing rates of deforestation, they have been ineffective in reforesting already degraded ecosystems as legal restrictions severely limit the use of forested areas (Legal Reserves and Areas of Permanent Protection) on all rural properties [47]. These regulations have been difficult to implement as reforesting lands is viewed negatively by landowners as the assumption is that once the forest recovers, the land becomes worthless or untouchable. The question remains as to how to incentivize transformative changes to the landscape when the predominant perception is that monoculture crops, with inevitable forest loss and detriments to the agroecological landscape, produce higher yields?

Aiming at reintegrating degraded agroecosystems into ecologically and economically functional areas, we developed and implemented Productive Agroforestry Restoration models focusing on restoring degraded or underused agricultural land into a multispecies productive system maintained as a forested environment [40]. The goal is to allow for a variety of outputs that take advantage of the inherent spatial and temporal variations of the system and produce direct (e.g., crops) and indirect (e.g., ecosystem services) benefits. We designed Productive Agroforestry Restoration as a response to the need for innovative productive systems that can generate income and restore ecosystems for the benefit of rural communities and society in general, offering land management solutions that can be implemented across the region.

Thus, the process of developing Productive Agroforestry Restoration models began with the premise that any land to be restored or (re)integrated into a sustainable agroecosystem (namely agroforestry) should not only focus on the restoration of ecological attributes but also be integral to the socioeconomic reality of the farm; the system must be both ecologically and economically sustainable. In 2011, with financial support from the Brazilian National Council for Scientific and Technological Development (CNPq) and EMBRAPA, we began implementing a project to leverage LEK in order to optimize traditional management practices and create agroforestry models for restoration. In collaboration with farmers, we developed a comprehensive list of potential agroforestry systems and species that farmers would be interested in cultivating. From this, we collectively chose models that were deemed more feasible and implemented these models at the EMBRAPA Research Station in Caçador (ERSC) replicated across an area of 40 ha. The results [35,36,39,40] allowed us to expand the implementation of such models to more than 50 small-scale farms in the region. Importantly, we also integrated the needs and expectations of farmers as four key requirements of the models. Specifically, the models must be easy to understand and implement; fast, through the rapid (re-)establishment of a forest canopy by using fast-growing pioneer species; profitable, as investment should result in economic return; and flexible to regional characteristics, property, and goals. Implementation can take place at different places and scales in a property and be integrated into a landscape restoration program.

The species selection for Productive Agroforestry Restoration varies depending on the region in which it is implemented. Initially, a few key species should be identified in order to fulfill the need for rapid forest cover establishment and acting as a cash crop. In many tropical and subtropical regions, the use of native, fast-growing species from the legume family is highly recommended for their multiple benefits, which include use for firewood and lumber, rapid deposition of soil litter, and nitrogen fixation, among others [59,60]. As some fast-growing trees are short-lived, sometimes with cycles of less than 20 years, their replacement should be considered early in the management planning. On the other hand, cash crops should be based on long-lived species in order to create financial predictability. There is a myriad of species and combinations but the design of a system with initial complex arrangements should consider possible constraints, such as seed/seedling availability, potential market outlets, local labor capacity due to systems with higher management requirements, and lack of technical knowledge about the species and their interactions, among others. Thus, biological complexity and product variety is often better achieved after the establishment of an initial simple system.

Among the Productive Agroforestry Restoration systems implemented at the ERSC and currently being used in farms across Southern Brazil, the most successful was designed to be implemented in any region of the subtropical highlands of Southern Brazil and is based on two species: the fast-growing bracatinga (*Mimosa scabrella*) and erva-mate. While bracatinga is a legume pioneer tree that regenerates spontaneously in the region and is used for firewood and slab props, erva-mate is a shade tolerant low maintenance tree with a consolidated market.

This system was first used to restore a degraded agroecosystem that after several decades of high intensity commodity monoculture had extremely impoverished and compacted soil. Following the Productive Agroforestry Restoration model, we planted the fast-growing bracatinga in rows 6 m apart (1.5 m distance within rows; Figure 4a) that rapidly formed a forest canopy, as shown in Figure 4b. As a pioneer species, bracatinga is expected to very rapidly form a forest canopy as it reaches average heights of 2.8 m after one year, 6.8m by year two, and 9.5 by year three (at which point a forest environment is established), and 13.3 m by year four (diameter at breast height of 2.2, 7.6, 11, 13.3 cm, respectively) [36,61,62]. With ideal shaded conditions for erva-mate in place with the establishment of forest cover by year three, erva-mate seedlings are planted in double rows between bracatingas (~3000 plants per hectare), as shown in Figure 4c.

**Figure 4.** Productive Agroforestry Restoration carried out in a degraded agroecosystem: (**a**) initial planting of the fast-growing bracatinga for rapid canopy development (year 1); (**b**) developed canopy allowed for the plantation of erva-mate (year 3); (**c**) commercial maturity of erva-mate generates profitability for the system (year 8); (**d**) different methods of bracatinga harvesting: left—alternate full row removal for optimized financial return from lumber; right—in row-alternate tree girdling (trees with brown crowns) for subsequent harvesting to minimize environmental change.

Bracatinga demands yearly pruning until at least year three to obtain trees with higher market value (timber without knots) and high branching that does not disturb erva-mate cultivation [40], while erva-mate should be cultivated following the silvicultural practices developed described above. Combined with the planting of bracatinga, in the first two years we cultivated soy between rows which helped to ameliorate revenue/investment financial ratio during this period. Later, by year five, income was again generated through bracatinga thinning (at 50%) and finally at year eight, erva-mate harvesting reached commercial levels when the system become profitable, as shown in Figure 4c, with cost–benefit ratio varying between 2.0 and 2.6 with higher levels obtained when bracatinga is harvested for lumber [36,40].

The system presented offers an alternative for land restoration, but its elements are open to variation and diversification. As environmental conditions gradually improve over time, especially soil structure and fertility, with much lower levels of humidity and temperature fluctuations, other species can be integrated to the model, taking advantage of the horizontal and vertical space available, which includes vines, herbs, and shrubs for various uses (food, medicinal, handcraft, etc.). Importantly, landscape managers can support the use of these systems as a means to reforest Legal Reserves on rural properties, while remaining productive.

We are currently undergoing a comprehensive monitoring and modeling effort in order to quantify and qualify the socioeconomic benefits and improvement of environmental services provided by the practices discussed herein. Our initial results show extensive improvement in terms of forest species diversity, which increased from the initially two tree species planted to 40 species successfully recruited [39] in five years. Furthermore, the occurrence of vines, shrubs, herbs, and an abundant natural regeneration in a multi-strata forest seems to confirm the widely held assumption that species diversity is expected to increase with habitat complexity [63]. Moreover, those results contrast sharply with our control area (no intervention since 2010; unpublished data) in which no trees to date have managed to recruit. In terms of the land sparing approach, our results suggest that merely setting aside land for forest restoration and conservation may be insufficient. Our model does not have some of the common pitfalls of other restoration systems that focus solely on maximizing substrate stability or primary productivity, often resulting in arrested succession and demanding additional efforts to encourage successional change [64]. While detailed carbon monitoring is underway, we have already been able to estimate large-scale benefits of restoring riparian forests in the region [65], which can be used as a proxy to estimate the benefits of restoring degraded lands and forests.

#### **4. Stakeholder Engagement and Current Challenges**

The introduction of new technologies in traditional agroforestry systems help initiate a process of social, economic, and environmental transformation in family farming involved in the participatory research. Firstly, our focus on participatory research that values the knowledge of small-scale farming families has led to self-reflection and a growing self-awareness of the value of this knowledge and the associated environmental identities in wider socioenvironmental discourses, and also their rights as farmers and food producers. Our focus on documenting and sharing knowledge across institutional, class, and gendered divisions has enabled the creation of a knowledge network that has led to many new initiatives, ideas, and solidarities. We are in the process, for example, of supporting the development of a network of women farmers that will identify the needs and challenges women face in participating as active members of decision-making and knowledge sharing circles, not only on the farm, but in the wider context of agroecological production. Farm visits that bring together a variety of farmers, practitioners, and researchers have also shown to be fruitful in knowledge exchange across various spheres, leading to innovative management strategies for pruning and pest management, among many others. Our initiatives are also helping to consolidate our partner communities into a collective network with a stronger political voice. One major advancement has been the creation of a strategic council (*Observatório dos Sistemas Tradicionais e Agroecologicos da Erva-mate do Paraná*) for traditional and agroecological erva-mate production systems, spearheaded by the Public Prosecutor of Labor of

Paraná, which supports the continuation and expansion of these systems. This council brings together 27 organizations that are working together to bring greater awareness to the ecological and cultural value of traditional erva-mate production, while also incentivizing new products, markets, and other economic benefits.

In terms of economic impact, optimization strategies, such as those tested and implemented in the *caíva* systems, are showing promise in terms of improved incomes for farmers [20], while the models being tested at the ERSC have shown possibilities of economic returns in relatively short periods of time [40]. The restoration model is also being rolled out through partnerships with industry to better test these productive systems at a larger scale and gain more concrete insights into the economic capacities of these models. Part of our ongoing research is to determine the indirect environmental benefits of these agroforestry systems, particularly in terms of carbon capture, water quality, soil health, and biodiversity, as support for the development of payment for ecosystem services models. As noted above, traditional erva-mate agroforestry systems have been shown to have significant levels of tree biodiversity [18,40], and despite consistent cattle grazing for several generations, the tree regeneration potential of *caíva* systems has been shown to be quite strong, with significant levels of diversity in terms of regeneration in comparison to those found in the Santa Catarina Forest Inventory [44]. Clearly, these productive forest systems have substantial environmental resilience and offer compelling strategies that can be implemented across the region. Preliminary results on ecosystem services, as noted above, are also showing greater potential for total carbon and nitrogen capture in erva-mate agroforestry systems than in monoculture erva-mate production areas. Nevertheless, more detailed data is required to continue to inform policy and regulatory frameworks.

Despite the benefits, several challenges still face the continuation and expansion of traditional and agroecological systems in the region. Current legislation related to forest management, for example, severely restricts silvicultural practices on private properties, with relatively arbitrary quotas placed on the number of trees that can be removed from native forests, while regulations related to livestock grazing and production in silvopastoral systems remain unclear. Although current legislation has been important in stemming the devastating loss of forests that occurred throughout the twentieth century, small-scale farmers feel disproportionately affected by the regulations, which has led to mistrust on both sides of the issue [38]. The oral history interviews conducted as part of our research have clearly underscored how tensions between small-scale farmers and environmental agencies have led many farmers to question the continuation of these systems as the current impasse seems insurmountable [33]. Yet, through research and advocacy in collaboration with farmers, changes are taking place, with environmental agencies such as Paraná State Environmental Institute (IAP) and the federal Brazilian Institute of the Environment and Natural Resources (IBAMA) participating in recent events and outreach activities organized as part of this project, and in the *Observatório*. Promisingly, these agencies are looking to update regulations and change legal restrictions based on the current state of forests in Brazil, supported by the data and experiences projects such as ours are sharing.

Although changes are taking place within policy circles, one of the biggest challenges we face is the inherent bias not only against forests, as they continue to be viewed as useless, which is directly related to current legislation that prohibits use, but also the culture of agricultural research and outreach agencies that are very much focused on conventional agriculture and mechanization and/or modernization at the expense of traditional knowledge and practices. Agroecological or traditional farming practices are generally excluded from agronomy courses in universities and technical colleges, and as such, the majority of outreach workers have little experience engaging with these alternative approaches. Despite the myriad policies that have been enacted to support small-scale family farming and agroecology/organic production in Brazil (see [66,67]), these policies have not necessarily trickled down to have clear impacts on small-scale farming communities, while others have reinforced the dominant model of intensification based on monoculture commodity crops. Nevertheless, recent developments occurring through ongoing engagement with environmental

agencies and other institutions are starting to show promising shifts in the top-down approaches to governance and agricultural outreach.

Despite national policy initiatives implemented in the last 10 years that have focused on drastically increasing the amount of data related to land use in Brazil (for example the Environmental Rural Registry (CAR) [47] and the National Forest Inventory [68]), federal and state land management planning is still in its infancy. One of the major challenges land management attempts face, however, stems from conflicting perceptions of productive land use. On one side, agriculture research and outreach agencies and large agrobusiness are pushing to modernize production through intensive, high-input monocultures, such as soy and corn. While these systems are seen as more productive, offering greater yields and thus higher income than traditional systems, there is a range of negative consequences including human health and food security, loss of farmer autonomy and increased debt, deforestation and loss of biodiversity and LEK, and impacts on water and nutrient cycles. On the other hand, environmental policies focus on protecting the remaining old growth forests and attempt to increase forest cover through restrictive laws that prohibit the use of forest resources, for example through Legal Reserves in which 20% of the property must be forested [47]. Furthermore, there remains an underlying assumption across both agriculture and environmental policy areas that native forests are not productive land, meaning either that the land must be deforested to become 'productive', or the burden of maintaining the land as forest (i.e., untouchable) falls on the landowner. This conflict is disproportionately felt by small-scale producers who use traditional agroforestry systems because it is exactly these farms that continue to maintain forest cover, which often extends well beyond the required 20%. Yet the law prohibits most types of forest management and agencies that monitor and inspect farms tend to administer fines with the onus on the landowner to prove they are within their legal limits, an unrealistic requirement for most small-scale farmers (for a full discussion, see [30]).

Our research is demonstrating that there is middle ground between these two competing land use perceptions that offer economic, cultural, and environmental benefits at both the local and regional scales. Implementing such models and practical approaches at the landscape scale can have dramatic impacts on the amount of land under forest cover, with the inherent returns of improved ecosystem services, biodiversity, and carbon sequestration. It can also bring about significant changes to the economic, cultural, and social value associated with forests and the products derived from them. By supporting and disseminating the knowledge, environmental subjectivities, and intangible heritage associated with traditional agroforestry systems, small-scale producers have an opportunity to create and capture niche markets that value ecologically and culturally significant products and processes. Thus, agroforestry systems using native species can be productive and economically, culturally, and environmentally viable, and through a transdisciplinary approach, land managers can work with traditional producers to develop best practices that can be implemented on farms across the region to have a real impact on sustainable land use on a larger scale.

#### **5. Conclusions**

The research presented herein demonstrates innovative approaches to documenting, valuing, and leveraging traditional agroforestry practices as a means to support diverse, resilient agroecological landscapes. While the focus of our research is on small-scale farms, the models we are testing show potential for scaling-out, offering promising alternatives that landscape managers can use to support sustainable land use and land cover change. The collaborative approach to research has been fundamental in this project as we are integrating and valuing different perceptions of agroecosystems and ensuring that communities are central to all aspects of the work, from defining research questions, to developing monitoring systems, and disseminating the results. Grassroots initiatives and locally adapted agroecological practices, such as those used in erva-mate agroforestry, are often ignored by research and government agencies, resulting in serious disconnect between overarching policy frameworks, such as the United Nations Sustainable Development Goals, and actual strategies that have the potential for transformational change at the landscape scale. Our work is attempting to bridge

this divide by leveraging the knowledge and practices small-scale farmers have been developing for generations before they are lost to the dominant paradigm of conventional agriculture.

**Author Contributions:** The authors contributed to the development, implementation, analysis, and writing of this study as follows: conceptualization and methodology, A.E.B.L. and E.R.N.; formal analysis, A.L.H., A.E.B.L., and E.R.N; resources and funding, A.E.B.L. and E.R.N.; writing—original draft preparation, A.E.B.L. and E.R.N; writing—review and editing, E.R.N., A.L.H., and A.E.B.L.; project administration, A.E.B.L. and E.R.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by EMBRAPA, grant number (16.16.05.002.00.00), EPAGRI, Programa SC Rural of the Government of the State of Santa Catarina, Social Sciences and Humanities Research Council of Canada, FLEdGE Partnership Grant, and Wilfrid Laurier University's Centre for Sustainable Food Systems. ERN was supported PNPD/CAPES.

**Acknowledgments:** The authors are sincerely grateful to the farmers and their families who have shared their stories with us. Particularly we thank Maria Izabel Radomski, Ednilson Pereira Gomes, Demerval Pessin de Farias, Bernardo Vergopolen, Carlos Urio, Arnaldo Soares, João F.M.M. Nogueira, Ricardo Gomes Luiz, Alessandra I. de Carvalho, and Robson Laverdi.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Pollination Potential in Portugal: Leveraging an Ecosystem Service for Sustainable Agricultural Productivity**

**Caroline Wentling, Felipe S. Campos, João David and Pedro Cabral \***

Campus de Campolide, NOVA Information Management School (NOVA IMS), Universidade Nova de Lisboa, 1070-312 Lisbon, Portugal; cwentling@novaims.unl.pt (C.W.); fcampos@novaims.unl.pt (F.S.C.); jcdavid@novaims.unl.pt (J.D.)

**\*** Correspondence: pcabral@novaims.unl.pt

**Abstract:** As urbanization and agriculture increase worldwide, habitats and food sources for wild pollinators are often fragmented or destroyed. As wild pollinators contribute both resilience and variety to agricultural fields, it is desirable to implement land management practices that preserve their well-being and ability to contribute to food production systems. This study evaluates continental Portugal for its change in suitability to host bee's pollinator species (*Apis mellifera*) from 1990 to 2018. It uses the InVEST crop pollination modeling tool and CORINE Land Cover, as well as parameterization to produce pollinator abundance and supply maps. These are generalized to municipality boundaries to provide actionable insights to farmers and policymakers and strengthen land management practices. It finds that the potential for pollination services is growing, with averages of both pollinator abundance and supply indices improving by 8.76% across the continental territory in 28 years. The study results are validated using another pollination index derived from a study that is based on expert opinion and field sampling in a sub-region of Portugal. This method of aggregation of model results and comparison of the percent difference by administrative boundary has the potential to better inform both policymakers and farmers about the pollination potential on a local level, as well as inspire interventions for future productivity.

**Keywords:** land use changes; wild bees; land management practices; validation; InVEST model

#### **1. Introduction**

Ecosystem services are natural processes from which human benefit, whether directly or indirectly [1]. Because natural "capital" (i.e., trees, atmosphere, carbon, information, nourishment, etc.) and human reliance on it is difficult to quantify economically [2] its value is often discounted in policy development. However, these services have tremendous effects on our wellbeing, resilience, and markets [3], making them a valuable addition to discussions about sustainable land management practices [4].

Pollination is one of these services from which humans reap significant benefit [5]. Though wild bees provide essential pollination services to both wild plants and crops alike [6], "agricultural intensification jeopardizes wild bee communities and their stabilizing effect on pollination services at the landscape scale" [7,8]. This type of loss has impacts on national economies. The estimated annual value of ecosystem services provided by wild insects and other animal pollinators (including pollination, dung burial, pest control, and wildlife nutrition) equates to more than USD 57 billion [9]. Other estimates project losses of USD 1.4 billion of the gross domestic product (GDP) between 2011 and 2050 in the US alone due to pollinator loss [4]. Insect pollination accounted for 35% of global food production in 2004 as well as 75% of crop types, [8,10]. Losses in crop pollinators are expected to affect the world supply of fruits, vegetables, oilseeds, and cotton, leading to direct and indirect effects on global commodity supplies and prices [4]. Worldwide declines of pollinators can catalyze similar trends in wild plant species [7]. These implications on both human well-being and environmental vibrancy necessitate the utilization of

**Citation:** Wentling, C.; Campos, F.S.; David, J.; Cabral, P. Pollination Potential in Portugal: Leveraging an Ecosystem Service for Sustainable Agricultural Productivity. *Land* **2021**, *10*, 431. https://doi.org/10.3390/ land10040431

Academic Editor: Diane Pearson

Received: 12 March 2021 Accepted: 15 April 2021 Published: 17 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

models that can characterize the effects of current trends, as well as predict and evaluate potential future scenarios [11].

Though there is a multitude of species responsible for the pollination of human consumable crops [7], many farmers employ domestic bees for managed pollination. However, the utilization of wild bees increases temporal stability as well as additional efficiency for certain crop species [6]. Though not strictly necessary, some vegetable species yield higher quality and more pest-resilient crops, in addition to improved seed production [10] after being visited by wild pollinators. Heterogeneous and organic fields are usually more suitable, both in terms of habitat appeal as well as in food resources, which may attract pollinators within their foraging ranges [4,12,13]. Understanding these types of behavior and interdependencies can improve the way farmers and policymakers adjust their practices to improve yields, as well as maintain sustainable supplies of pollination services into the future [7].

Previous studies on pollinator suitability have been performed in sub-national regions around the world [13] at more general continental or global levels [2], are limited to specific crop types [14], or landscape types [10]. Some of the previous literature provide frameworks for incorporation into future studies [5], some leverage or derive theoretical monetization models [15] and some do not incorporate spatial dependence into their models [9]. Derivation and validation of these studies range from labor intensive field sampling [16] or leveraging of primary sources [7], to expert opinion [17], to predictions of models derived from environmental inputs (such as Land Use Land Cover (LULC), climate, or topology [8]).

This study demonstrates the viability of applying a spatially dependent model of pollinator suitability to an entire continental area (corresponding to the mainland area of Portugal, designated by the Nomenclature of Territorial Units for Statistics (NUTS) level 1 code PT1: "Continente" in Portuguese, or "continental" in English) and then aggregating the interim results to subregions (NUTS subdivision 3) to evaluate both overall trends and local changes over time. This aggregation provides new opportunities for land management and innovation practices applicable to various levels of local administration. This type of investigation has not yet been applied to Portugal.

To this end, the study evaluates the changes in pollination suitability in continental Portugal from 1990 to 2018 of a representative guild characterizing the behavior of the European honeybee (*Apis mellifera*). It derives pollinator abundance and supply indices from input LULC raster maps for the area as well as parameterized pollinator guilds. These are processed by the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) crop pollination model, which incorporates spatial dependency of nearby floral resources and nesting sites in relation to pollinator foraging ranges. The resulting raster maps are aggregated to administrative municipalities, and the percent variation (PV) is calculated. The results are evaluated for their trajectories of change and validated via the extrapolation of a local pollinator index based on in-field sample collected data and expert opinion.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Portugal is a European country of about 92,212 square kilometers on the southwestern corner of the Iberian Peninsula [18]. It contains the most western point in continental Europe and shares a land border with Spain. Portugal experiences Mediterranean climate of dry, hot summers and wet, cool winters [19], though this varies throughout the territory's microclimates (generally categorized as cooler and rainier in the north while drier and hotter in the south). According to CORINE Land Cover (CLC) of 2018, approximately 3.83% of the country's land cover is artificial surfaces, 47.81% is agricultural land, and 46.48% forests, with the remainder, made up of wetlands and water bodies (Figure 1) [20]. 28 years prior, artificial surfaces only covered 1.9% of the land, with 47.80% and 47.92% of the area dedicated to agricultural and forest land, respectively. This sizeable increase in artificial cover is consistent with the high rates of urbanization seen in and around Portugal's major cities. Desertification, the degradation of dryland also affects the changing classifications of land cover especially in the interior of the country [21].

**Figure 1.** CORINE Land Cover of 2018 (level 1) distribution in Portugal, with municipal boundary definition (Carta Administrativa Oficial de Portugal: CAOP 2018).

Portugal is composed of 308 "concelhos" (municipalities, NUTS 3), with 278 of these located on the mainland [18]. Of Portugal's 240.7 billion USD GDP in 2018, 2.05% was produced by the agriculture, forestry and fishing industry [22]. With almost half of continental Portugal's land surface devoted to agriculture, there is tremendous value in ensuring the successful production of cultivated crops. It is estimated that Portugal is home to more than one thousand pollinating insect species, including a variety of bee, hoverfly, butterfly, and flower beetle species [16,23,24]. As of 2018, 680 distinct bee species have been collected and catalogued in Portugal [25]. Within the River Minho area alone, 200 distinct species were catalogued for a smaller scale study on pollination services [16]. According to the Joint Research Centre (JRC) Technical Report, Portugal demonstrated the "highest increase of pollination potential" from 2000 to 2006 in the European Union (EU) [26]. The main pollination season in Portugal can range from March to September, which includes the season of most active airborne pollen particles in the country as well as the general period for crop pollinator foraging in the north half of the globe [19,26].

#### *2.2. Software and Data Management*

This study was carried out using the free and open source InVEST crop pollination model, available under the open data license [27,28]. It also leverages the proprietary ArcGIS software (ArcMap 10.6) to perform the spatial temporal variation model, visualize the results, and for validation. All data included in the study is open data freely available to the public through the portals described in Section 2.3.

#### *2.3. InVEST Crop Pollination Model*

InVEST is a software platform of the Natural Capital Project and a suite of models to evaluate and chart a variety of ecosystem services, ultimately to inform decisions on how to manage these natural resources by quantifying their economic impact. It has been used in previous academic studies to marry macroeconomic scales with local environmental processes to predict multiple future scenarios of varying degrees of environmental action to "resonate with political economy audiences" [4].

The InVEST crop pollination model produces pollinator abundance and supply indices, which are scaled from 0 (least suitable) to 1 (most suitable) [27]. Abundance index represents the likely location of their activity, while supply index describes the likelihood, based on proximal nesting sites and food resources of the location and foraging ranges of the species, for pollinators to nest in a space. The results characterize wild bee pollinator guilds (groups of bee pollinators demonstrating similar nesting and foraging preferences as well as foraging distances and relative abundance).

The model utilizes land cover raster maps as well as persistent bio and guild tables as inputs. It incorporates habitat parameters (estimated nesting site and floral resource availabilities, and relative abundance per guild) for each cell of the input raster, considering the floral parameters of its neighbors [20,26]. One of the key features of the model is the incorporation of foraging distance, which allows the model to bridge the possible spatial separations of nesting and foraging habitats [8,12,29]. This model was selected for its accommodation of the spatial dependency required in such geographically explicit studies.

To make meaningful comparisons between time frames, pollinator abundance and supply indices for each pixel were generalized via zonal statistics into the 278 municipalities under study. The resulting statistical means of each municipality were utilized to calculate the variation from 2018 to 1990, according to Equation (1):

$$PV\_{\varepsilon} = \frac{\Delta A bundance\_{2018c} - \Delta A bundance\_{1990c}}{\Delta A bundance\_{1990c}}\tag{1}$$

where *PVc* is the percentage variation index for delivering pollination abundance for year 2018 in comparison to the baseline year 1990 for each concelho©. The general flow is depicted in Figure 2, with a comprehensive modeling workflow for the percent variation of abundance index. This process was executed in ArcGIS software (ArcMap 10.6).

**Figure 2.** Schematic representation for analyzing the statistical comparison model flow.

#### *2.4. Spatial Data*

Land cover raster maps utilized are available from the Copernicus project, provided by the European Environmental Agency (EEA) [20]. The CLC classification includes 44 distinct subcategories that fall within five major areas: artificial surfaces, agricultural areas, forest and semi-natural areas, wetlands, and water bodies. As this study seeks to understand the change of ecosystem services over time, the earliest and latest available years (1990, 2018) are used. The raster data have a spatial resolution of 100 m and a minimum mapping unit of 25 ha [20]. All "slivers", landmasses associated with continental Portugal but removed by water, have been excluded from the study area. All data in the study is in the common coordinate system ETRS\_1989\_Portugal\_TM06.

The biophysical table (required for InVEST crop pollination) corresponds to the LULC classifications to establish suitability for nesting and floral resources of each raster input (see [27] for additional details). Of particularly high suitability are certain agricultural areas and forest edges, both of which tend to provide heterogeneity of habitat in the form of diverse nesting space and floral resources within a small area, often promoting insect activity [8]. Nesting and floral resources parameter values are provided in the supplementary material of [8]. The values are derived from expert opinion and leveraged in their European continent level of pollination, also utilizing CLC input raster maps.

This study utilizes the CLC classification conversion provided in a study on pollination services across the European continent [8], as it is directly applicable to the study area (Appendix A). The conversion parameterizes all 44 classifications of CLC, generalized as a single season (versus representation of seasonal pattern variations spanning a calendar year) and a single nesting substrate (no distinction between cavity or ground preference), thus the results are representative of these generalizations.

Guild parameters assign values to represent the different behavior patterns of various bee species. These patterns include preferences for different nesting sites and floral resources, as well as relative prevalence and foraging ranges (Appendix A).

Though InVEST has the capability to model multiple nesting types, seasons, and bee guilds, insufficient data exists in previous literature to leverage the full potential of the tool (let alone the variation of parameter values due to environmental conditions [12]). Therefore, values for individual parameters were aggregated from a variety of sources [8,14,30] to describe one pollinator guild. *Apis mellifera*, better known as the European honeybee, is considered "the most economically valuable [pollinator] of crop monocultures worldwide" [8] and is widely employed in managed crop pollination and honey production and is native to mainland Portugal. This species is well studied and can be easily characterized as per requirements of the InVEST model.

The Carta Administrative Official de Portugal 2018 (CAOP 2018) is originally available from DGT (Direção-Geral do Território), a portal providing geodesic and geographic information services by the Portuguese ministry of agriculture, sea, and environment, and territorial management [2,11,13]. It includes 278 continental Portugal administrative territories, ranging in size from 7.94 to 1720 km2 (São João da Madeira and Odemira municipalities, respectively).

#### *2.5. Validation*

The InVEST ecosystem service modeling toolset is well established and widely used in academic study [15,27], which supports its reliability. However, it has some recognized shortcomings and is subject to the quality of input data [16]. Validation of the results is required prior to their influence on future decisions on the management of ecosystem services.

So that the study results can be meaningfully compared to the validation methods, the pollination indices are normalized such that they are distributed between their reported minimum index (adjusted to 0) and their reported maximum value (normalized to 1). Another study performed in a subsection of continental Portugal is leveraged as validation. The study developed a Pollination Suitability Index for Riverine Landscapes (PSIRL) in the River Minho (norther border of Portugal with Spain) in 2018 [27]. Though this approach has

its own limitations (the study considers insect pollinators in general, not just *Apis mellifera* and it derives specifically for riparian areas, requiring generalization and translation to the input CLC LULC), its index is derived from expert judgment, floral diversity, and actual field surveys, increasing the overall confidence in the results. No other spatially comprehensive yet reliable data exist in Portugal for validation.

The validation process is depicted in Figure 3. The original PSIRL index is translated to land use codes used by the CLC LULC (see Appendix A), The "unclassified" features representing water areas are selected then intersected to create water edge lines. These are buffered by 10 m and merged with the LULC polygons. These are then converted back to rasters and zonal statistics are applied. These are joined to the CAOP municipal boundaries, normalized, and then compared to the normalized index values of the results.

**Figure 3.** Pollination Suitability Index for Riverine Landscapes (PSIRL) validation flow.

#### **3. Results**

#### *3.1. Land Use Land Cover Evolution from 1990 to 2018*

Table 1 demonstrates the land surface utilization (LULC areal data) within Portugal in 1990, and then the percent variation from this baseline to 2018. These have been generalized to the broadest category (CORINE Land Cover designation level 1). Note that each L1 category is not associated with homogeneous biophysical parameters. The table demonstrates a doubling of artificial surfaces between 2018 and 1990 (largely inhospitable to bees), as well as a 1.8% and 3.0% drop in largely appealing habitats (agricultural and forest cover, respectively) over time. This would suggest an overall decrease in pollination services over time.

**Table 1.** Land Use Land Cover (LULC) percentage by 1990 and its percent variation (PV) by 2018.


#### *3.2. Pollinator Abundance and Supply Indices*

Figure 4 depicts the spatial variation of pollinator abundance and supply indicators for 2018. As one would expect, both indices follow similar spatial patterns: less hospitable in the urban and water areas (red), more hospitable in forested areas (green). Coastal areas are particularly unfriendly, both in the West and the South, with much of eastern Alentejo region exhibiting low suitability as well. On the other hand, much of north eastern Portugal and around the border of Alentejo and the Algarve regions appear to be quite suitable along both indices.

**Figure 4.** Spatial distribution of InVEST crop pollination model outputs (**a**) abundance index and (**b**) supply index in 2018.

Table 2 displays the general statistics of the raster results. The minimum index value remains zero for those areas that are unsuitable for pollinator activity (both inhospitable to bee nests as well as outside of the range of foraging). Overall, the means and maximum values for each index have increased between 1990 and 2018, indicating that the overall suitability for pollination services in Portugal is growing. This yields a slightly larger standard deviation, which reflects a greater range of index values distributed across continental Portugal.

**Table 2.** Statistics of InVEST crop pollination model results in 1990 and 2018.


*3.3. Pollination Service Changes from 1990 to 2018*

Once the results are associated by municipality, inferences about trends for each administrative boundary are more easily understood. Ideally, this will contribute to better policy making at the district level. Figure 5 displays both the percent variation of abundance

(color of polygon area) and supply (color and size of overlaying triangle,) from a baseline of 1990 to 2018. Red colors indicate negative changes in suitability indices, while green and blue indicate positive trends. Yellow indicates no significant change over the 28-year study period. The triangle size also indicates the degree of deviation of supply from the 1990 baseline. Both changes in supply and abundance tend to fall into the same categorizations. Those areas that experience differences are usually (but not exclusively) characterized by a supply index that is slightly more extreme than that of abundance. This suggests that the pollinators may be more selective about their habitats (origins) than in the areas they are willing to traverse in search of food.

**Figure 5.** Percent variation (PV) from 1990 to 2018 of abundance and supply indices by municipality.

#### *3.4. Validation*

The PSIRL validation technique required translation of the given PSIRL index values to the original input LULC raster map (Appendix A). The resulting difference map (Figure 6) includes an overlay of the riparian areas (including a buffer of 300 m from water areas), as well as an indication of the area of which the PSIRL index was originally derived (the River Minho area in north western Portugal, identified with a red box). The yellow zones indicate those in which the validation demonstrates good coincidence with the results of the study (within a 5% tolerance), whereas the stronger purple and red colors indicate larger discrepancies between the two methods (a maximum discrepancy of 38%). Clearly the results incorporate some amount of spatial autocorrelation, though surprisingly these are not necessarily correlated with riparian adjacent areas as one might expect. The PSIRL tends to slightly over-predict the supply of pollinators as compared to the results of the study. The difference map indicates a strong spatial similarity between the two models, strengthening the confidence in the study results.

**Figure 6.** Difference between study results and Pollination Suitability Index for Riverine Landscapes (PSIRL) with riparian zone overlay. Red box indicates River Minho area, in north western Portugal.

#### **4. Discussion**

#### *4.1. Study Significance*

Bees require suitable places to nest and sufficient food sources near nesting sites to sustain them [1]. These and other factors have been applied to a model that produces maps of projected pollinator activity within Portugal. Too often, stakeholders (farmers, policy makers, economists, etc.) ignore the subtle interactions between ecosystem services and production, which can be to their own detriment when those nebulous costs outstrip their values [4]. It is estimated that, at the current rate of land use transformation, the United States gross domestic product (GDP) will suffer a loss of 0.02% (or 15 billion USD) due to reduced wild pollinator habitats near agricultural sites [8], which can have ripple effects in other industries to compensate for the deficit.

Results of this study and other such investigation will ideally support farmers, land developers, and policy makers alike with better information from which to make decisions about how to better manage these resources as well as improve economics systems that depend on them by maximizing their sustainability. Agriculture and thriving pollinator communities are not mutually exclusive. In fact, well-managed cropland can be economically and ecologically productive [6,7,10,11]. For instance, farmers could identify locations for crops based on maximizing exposure to wild pollinators, adjust their management towards organic practices, or maintain heterogeneous nesting substrates that would attract diverse and productive pollinator populations [13]. Likewise, configuring farms towards a variety of pollinators (instead of just the domesticated varieties) can produce better yields, as different pollinators are associated with varying levels of productivity for certain crops [29]. Even the understanding of the tendency of larger bees to populate new fields and smaller bees to prefer older fields [8], or the observed abundance and variety of pollinators in forest edges and grasslands [6] can assist with the development of management strategies. Further, the understanding of the relationships between space, crops, and pollinators may

provide an incentive to better care for areas beyond crop fields [31] or mitigate the appeal of monoculture practices [7]. In fact, there is potential for the benefits of ecosystem services management practices to positively impact other areas within foraging distances of the appealing habitat sites. On the flip side, poor planning regarding conversion patterns of forest to agriculture can have devastating impacts on wild bee populations that will also undercut the productivity of the new agricultural land [7].

The association to municipalities provides a simplification of the detailed information to ease comprehension of the big picture, such that areas requiring intervention (those tending towards lower suitability) can be triaged and evaluated more efficiently. Though pollination may not be strictly required to achieve sufficient caloric intake, indeed many staple foods do not require this type of sexual reproduction, the production of many valuable nutrients require pollination [27], and pollination services haves been linked to qualitative (nutritional content, appearance) and quantitative (production yields) factors that boost economic value of agricultural production [32]. Ideally, stakeholders will be able to model and evaluate different policies and their effects on farm productivity, optimizing both resilient biodiversity as well as economic yields [1,6,10]. Methods to achieve this could be coordinating reserved land areas that provide pollination services through integration of natural areas throughout agricultural areas [31]. Further, incentive programs that promote healthy management practices or payment schemes could be organized, in addition to the inherent benefits experienced by the implementing farmers [26]. As this is a relatively new concept, there is much room for novel methods of accounting for ecosystem services within the economic structure of farmers and other land managers.

#### *4.2. Critical Analysis*

This study provides a valuable baseline indicator of pollinator services within continental Portugal. Overall, since 1990 there have been significant, polarizing tendencies of municipalities across Portugal. The percent variation of likelihood for pollinators to be active between 1990 and 2018 swings from −69% (Pedrógão Grande) to 107% (Ponte da Barca and Vila de Rei). As one would expect, there are concentrations of negative trend areas that are associated with major city areas and likely rapid urbanization (such as Lisbon, and Porto areas), as well as some areas of vast agriculture swaths (the south west portion of Alentejo Central), which are less hospitable. On the other hand, it is promising to see the constant improvement through central northern Portugal. Interestingly the areas of greatest percent increase and decrease are located adjacent to each other: the municipalities along the border of Médio Tejo and Beira Baixa both exhibit extremely positive trends since 1990, yet just across the border, several municipalities in Região de Leiria include some of the most negative changes in the same time frame. This is due to recent fires resulting in large swatch of burnt areas in these regions, making them inhospitable to pollinators, though their neighboring forested areas demonstrate favorable habitats. Some areas have significantly changed from the 1990 baseline. Ponte da Barca continues to improve its tendency towards pollinator likelihood in both abundance and supply, rising to a high of 107 and 109 percent variation increases, respectively.

The positive trend of pollinator abundance and supply indicators are consistent with the findings of another study that Portugal demonstrated impressive improvements in pollination potential [7]. From the maps of Figure 5, policy makers and farmers alike may better understand the existing trends of pollinator suitability since 1990, using this information to support new interventions that may increase pollinator suitability within each municipality. Of course, pollinator suitability may not be a priority to certain urban areas such as Lisbon or Porto, or municipalities specifically cultivating crops that do not require pollination such as the Douro region. Other areas of agricultural swatches that demonstrate reduced or unchanging suitability may benefit from a re-organization of agricultural land to better suit natural pollinator activity. These include areas such as Alentejo Central and Algarve areas, though many other persistent or worsening regions are distributed throughout Portugal.

#### *4.3. Additional Findings*

The results also demonstrate the stark impacts of forest fires on pollinator suitability, such as the dramatic changes in the Região de Leiria. Portugal is prone to fires in the hot, dry summers. Though these are often uncontrollable natural phenomena, understanding their effects on pollinator activity among other ecosystem services in conjunction with social loss may strengthen the attempts to better manage forest areas and inspire more radical interventions to recover the areas in the wake of such devastation. More granular studies may consider excluding burnt areas from their studies, though they were retained here as they contribute to the overall trends (encompassing both natural and human influence) in mainland Portugal.

#### *4.4. Research Limitations*

Though the results of this study are promising, there are several limitations of note and opportunities for future improvement. The results of this study are limited the to the available data and certainly leave much opportunity for further evolution. Because pollinators can differ significantly from ecosystem to ecosystem [16], leveraging the parameters of similar studies in other regions is often inappropriate. Better characterization of local bee species throughout the study area may yield more accurate characterization of the potential of this ecosystem service. For example: the pollination potential characterized in this study is relative only to *Apis mellifera*, which has different habitat and foraging preferences and activities (such as potential foraging distances) than other smaller, wild species. However, due to scarcity of data on the behaviors and preferences of other wild bee guilds in Portugal, only a single bee specie was characterized. The application of the methodology undertaken by [26] in the Minho river area (counting the number and characterization of pollinators active in a particular area) to the entire country was outside of the scope of this project but this and the inclusion of expert based models (EBM) could enhance future research [27]. In addition to better characterization of pollinators, more detailed parameterization of the biophysical table of LULC designations (such as the inclusion of multiple nesting substrates and seasons) could more accurately reflect the actual pollination activity throughout the year.

InVEST models measure the potential of the study area to provide pollination for bee pollinators. Additional considerations outside of the model purview will affect the actual pollination supply, such as the lack of accounting for pollinator persistence over time. Likewise, many other non-bee pollinators (such as butterflies, bats, moths, and birds) that are active in Portugal are not accommodated in the model. These and other such inherent limitations are described in the InVEST documentation in greater detail [11]. Notably, the model does not distinguish between natural or artificially initiated changes in pollination potential-discerning the source requires savvy technicians and good understanding of the local context to presume.

Regarding the input raster maps, the minimum mapping unit of 25 hectares of the CLC LULC data does not accommodate the impact of potential pollinator habitats or foraging supplies smaller than this area (such as in green spaces in urban areas). Similarly, the study does not accommodate the implementation of agricultural practices that may alter the desirability of the area for pollinators, such as the accommodation of nesting sites or use of pesticides. Though pollination is sensitive to both aggregation and spatial resolution as an ecological service that involves stocks and dynamics, it is expected that the CLC mapping units are appropriate for the scale of study [27]. Likewise, the study is subject to the accuracy of the CLC classifications. Any assumptions or misclassifications will propagate through this study.

#### *4.5. Future Opportunities*

The InVEST crop pollination model, in additional to providing suitability indices of pollinator habitat and foraging supply, can model a yield index for pollination impacts on existing agriculture [8]. This requires vector data detailing the geospatial location of farms along with their crop types, dependence on pollinators, abundance of managed pollinators, farm nesting sites and floral resources. This study did not have access to national farm data, and the establishment of the required attributes would require additional investigation outside of the scope of this project, likely including the need of expert opinion to properly assign values to these parameters. This is an area of potential study in future endeavors.

Suitability of edge environments has been noted to differ from that of non-edge environments. For example, forest edges are particularly suitable as pollinator habitat [13,16], but are not accommodated in the InVEST model. Other studies have included additional characterization of these areas, which could improve the approach with this additional nuance. Roadsides and riparian areas are further examples of opportunities for model adjustment and finer characterization. The study could also benefit from the accommodation of more granular parameterization, including that of urban areas hosting managed bee colonies or integrating green and biodiverse areas within the built environment, distinction between forest compositions, or refuge areas that may experience different pollinator assemblages.

#### **5. Conclusions**

The land use land cover impact on pollination services distributed over time and space was studied in the context of continental Portugal. The InVEST crop pollination service modeling tool was leveraged to understand the spatial relationship between pollinator abundance relative to a landscape's available habitat and food resources in accordance with their behavior and preferences. The results demonstrated an overall improvement in wild pollinator hospitality across the country, though several municipalities are becoming increasingly weaker in their suitability for such services. The relative distribution of pollinator hospitality indices was validated via a local pollination index based on field sampling and expert opinion.

In Portugal the measured distribution of tons of crop production in the country is almost equally distributed between known dependency (34.2%), non-dependency (33.4%) and unknown dependency on pollinators [32] With more than a third (and potentially up to two thirds) of the production weight relying on pollinators, there is a large economic incentive for agro-farmers, the primary beneficiary of pollinator services [33], to incorporate pollination ecosystem services into their practices and protect these resources. Further, secondary beneficiaries—such as consumers of more nutrient dense crops or governments receiving greater tax revenues—will experience positive effects from the products of these measures as well.

**Author Contributions:** Conceptualization, C.W., F.S.C., J.D. and P.C.; methodology, C.W., F.S.C., J.D. and P.C.; formal analysis, C.W.; investigation, C.W., F.S.C., J.D. and P.C.; writing—original draft preparation, C.W.; writing—review and editing, C.W., F.S.C. and P.C.; visualization, C.W. and F.S.C.; supervision, F.S.C. and P.C.; project administration, P.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported through the FCT (Fundação para a Ciência e a Tecnologia) under the projects PTDC/CTA-AMB/28438/2017—ASEBIO and UIDB/04152/2020—Centro de Investigação em Gestão de Informação (MagIC).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available on reasonable request from the corresponding author.

**Acknowledgments:** The authors would like to thank the anonymous Reviewers for their comments which helped improving the quality of this manuscript.

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

#### **Appendix A. Model Parameterization**

**Table A1.** Land Use Land Cover (LULC) classes, parameterization of the biophisical table for the InVEST model (nesting and floral resources), and PSIRL final scores used in validation.


#### **References**


### *Article* **Modelling the Impacts of Habitat Changes on the Population Density of Eurasian Skylark (***Alauda arvensis***) Based on Its Landscape Preferences**

**Nándor Csikós and Péter Szilassi \***

Department of Geoinformatics, Physical and Environmental Geography, University of Szeged, Egyetem utca 2-6, H-6722 Szeged, Hungary; csikos@geo.u-szeged.hu

**\*** Correspondence: toto@geo.u-szeged.hu

**Abstract:** The dramatic decline of the abundance of farmland bird species can be related to the level of land-use intensity or the land-cover heterogeneity of rural landscapes. Our study area in central Europe (Hungary) included 3049 skylark observation points and their 600 m buffer zones. We used a very detailed map (20 × 20 m minimum mapping unit), the Hungarian Ecosystem Basemap, as a land-cover dataset for the calculation of three landscape indices: mean patch size (MPS), mean fractal dimension (MFRACT), and Shannon diversity index (SDI) to describe the landscape structure of the study areas. Generalized linear models were used to analyze the effect of land-cover types and landscape patterns on the abundance of the Eurasian skylark (*Alauda arvensis*). According to our findings, the proportions of arable land, open sand steppes, closed grassland patches, and shape complexity and size characteristics of these land cover patches have a positive effect on skylark abundance, while the SDI was negatively associated with the skylark population. On the basis of the used statistical model, the abundance density (individuals/km\*) of skylarks could be estimated with 37.77% absolute percentage error and 2.12 mean absolute error. We predicted the skylark population density inside the Natura 2000 Special Protected Area of Hungary which is 0–6 individuals/km\* and 23746 ± 8968 skylarks. The results can be implemented for the landscape management of rural landscapes, and the method used are adaptable for the density estimation of other farmland bird species in rural landscapes. According to our findings, inside the protected areas should increase the proportion, the average size and shape complexity of arable land, salt steppes and meadows, and closed grassland land cover patches.

**Keywords:** land cover; land use; landscape structure; Eurasian skylark; farmland birds; prediction; Natura 2000

#### **1. Introduction**

In the terrestrial ecosystems of the world, the dominant land-cover category is agriculture (38%), including the arable-land use type [1]. In Europe, this value is much higher, at 45% (EBCC, 2015). The agricultural land-cover category contains various land-use types with different levels of human impact. The heterogeneity and spatial structure of these land-use/land-cover (LULC) patches vary greatly across rural areas, which has strong impact on farmland-bird diversity in Europe [2,3]. Many articles have determined that the decreasing trend of farmland birds is strongly connected with the intensity of agricultural management (level of use of fertilizers etc.) [4–7]. Very few studies have investigated the dramatic decline of the abundance of farmland birds, and its connection with change in landscape structure and land-cover heterogeneity [7–9]. There are some regional (country) scale studies that analyze the connection between land-cover types and farmland-bird population data [10–15]. These studies have indicated that the abundance of farmland birds is significantly connected with the intensity of agricultural cultivation, crop heterogeneity, and land-use change. Most articles focus on small, local study areas and analyzing the connection between Eurasian skylark (*Alauda arvensis*) abundance, and the proportions of crop

**Citation:** Csikós, N.; Szilassi, P. Modelling the Impacts of Habitat Changes on the Population Density of Eurasian Skylark (*Alauda arvensis*) Based on Its Landscape Preferences. *Land* **2021**, *10*, 306. https://doi.org/ 10.3390/land10030306

Academic Editor: Diane Pearson

Received: 15 February 2021 Accepted: 16 March 2021 Published: 17 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

type, height, coverage and heterogeneity [4,6,10,16–19]. The skylark does not prefer the fragmented landscapes by urbanized area, road network, hedgerows and heterogeneous land cultivation areas [7,20]. The agriculture is the dominant land use (matrix) of the European NATURA 2000 network, where the size and shape characteristics of different LULC patches, and the land cover heterogeneity can be essential for the protection of farmland bird species. Therefore, we hope that our results can be adding some new suggestions for the landscape planning and habitat design of national parks, NATURA 2000, and other protected areas. Our research also can provide important component for achieving the goals of the EU Birds directive [21].

The skylark is one of the most common farmland bird of rural landscapes in Eurasia, including Hungary. In the European Union, the Eurasian skylark has a declining trend in population between 2000 and 2018: Norway −47%, Lithuania −41%, France −38%, Czech Republic −29%, Hungary −24% and Germany −17%. Most individuals that breed in Central Europe spend the winter in the Mediterranean region, but small groups can stay in Hungary for winter [22]. This bird species have been introduced into the Nearctic, Australia and New Zealand [23,24]. From large-scale studies, habitat preferences, including for crop structure and heterogeneity are well-known. On the basis of small-scale regional-level studies, the regional-scale habitats and land-cover heterogeneity preference of a given species can be understood [10]. However, the connection between the spatial pattern of LULC patches (described with landscape indices), and skylark abundance is not clear.

In this study, we describe the landscape structure of rural landscapes with a very detailed (20 × 20 m minimal mapping unit) LULC map, the Hungarian Ecosystem Basemap (HEB). Comparing skylark abundance data with the HEB, we could identify preferred and nonpreferred skylark habitats, and calculate their landscape indices. The preferred habitat was separated into arable lands and grasslands because we wanted to analyze the effect of arable land and grassland landscape metrics on the skylark population. According to the pattern and process paradigm, which analyze the relationship between the landscape patterns spatial distribution and landscape processes, landscape indices are widely used as indicators of biodiversity and habitat changes [13,25–28]. After we identified preferred and non-preferred habitats for skylarks, we could calculate shape- and size-related class-level landscape metrics, and land-cover heterogeneity, and estimate the collective impact of these variables on skylark abundance [9–11,13,29,30].

The main goals of this study were to:


According to our hypotheses the population density of skylark is predictable based on the preferred LULC categories of skylark and landscape indices (proportion of LULC categories and shape and size related landscape metrics). The methodology is adaptable for analyzing the impact of landscape composition on other farmland-bird populations, and for predicting the population density of the skylark, in protected areas, where field observation-based datasets are not available.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Hungary is located in the Carpathian basin (45◦43 to 48◦35 N and 16◦06 to 22◦53 E) in central Europe, and is part of the Pannonian biogeographical region (Figure 1). The total area is 93,033 km\*, and its elevation ranges from 77 to 1014 m a.s.l. The most important land-cover type (61%) is agricultural land [31]. A further 20.7% is natural and seminatural grasslands and forest, and 5.5% is built-up area. In the 1990s, a dramatic landscape change was mainly caused by land privatization. Agricultural lands with low quality and poor

agroecological conditions were abandoned [32]. The common agricultural policy of the EU (strong decline of grazing livestock) and land abandonment caused the transformation of arable lands into non-cultivated lands, and the fast and spontaneous reforestation of open grasslands [33].

**Figure 1.** The spatial distribution of the MMM survey observation points in Hungary, where the Skylark occurred in 2015 (3049 observation points).

#### *2.2. Databases*

#### 2.2.1. Skylark-Abundance Data

In Hungary, a countrywide bird-monitoring survey has been conducted every year (like in 2015) by approximately 800 field surveyors who add their field-observation datasets into the Hungarian Common Bird Monitoring Database (MMM) [34–36]. The volunteers were not randomly distributed across Hungary. The survey allowed that the observers choose their area of observation. Each observation point received two spring visits, and the abundance of birds was observed (by hearing and visually) within a 100-m radius of each point. There is a minimum 500 m distance between the observation points. The surveyors left a minimum of two weeks between visits in mid-April and mid-June. The count was accomplished between 5:00 and 10:00, when wind speed was less than 5 m/s and there was no rain. Each observation point contains the average number of observed birds which were counted at the point in the two spring visits [34,35]. In 2015, surveyors counted 6763 skylark individuals across 3049 field observation points (mean value: 2.22, maximum: 34, standard deviation: 4.38.). We used MMM survey points from 2015 in the study area because the HUB land-cover map was also available from that time scale. We analyzed the proportion and spatial configuration of the landscape in the 600 m radius surrounding of the MMM observation points. 600 m buffer zone was chosen, because many author found that landscape composition and land cover types have the highest impact on the abundance of this species within this radius [10,37]. Land use types also have an effect on abundance of skylark population within 600 m buffer radius [38]. We used a very detailed (20 × 20 m mapping units) country scale LULC HEB maps [39] for analyses of the LULC characteristics inside these buffer zones. Unfortunately, the more detailed country scale statistical datasets about the crop structure surroundings of the observation

points were not available. Most of the MMM observation points (43%) is situated inside the NATURA 2000 SPA Protected areas, where the grasslands are mowed one-time every year after 15th of June.

#### 2.2.2. Land-Cover Database—Hungarian Ecosystem Basemap

The digital LULC HEB was created by the Hungarian Ministry of Agriculture. The basis year of this database is 2015. This very high resolution LULC dataset was based on other LULC maps of the European Copernicus Program, such as Urban Atlas, Corine Land Cover and High-Resolution Layers, and Sentinel-2 images. The dataset has a 20 × 20 m resolution (minimal mapping unit) and three category levels. Six classes in Level 1, 22 classes in Level 2, and 56 classes in Level 3 (see Appendix A Table A1). The database also contained three additional LULC categories in Level 4. We used the second level for analysis, and regrouped the LULC classes to reduce the number and the likelihood of autocorrelation between them. Our dataset for statistical analyses contained the following main LULC categories inside the buffer zones (Figure 2.):

**Figure 2.** Proportion of the main land cover categories in the 600-m buffer zones, where the Skylark abundance were detected (3049 observation points) based on Hungarian Ecosystem Basemap.

In our investigations, we aggregated the LULC categories of the HEB database, such as "forest", "wetlands and water surfaces" LULC categories (Table A1). The HEB web map and its documentation are freely available (downloadable) on this website: http: //alapterkep.termeszetem.hu/ (accessed on 15 February 2021). [39]. 46% of the country is arable land and cereals take the 62% of the arable lands. According to the country scale statistical datasets, the proportion of the crop structure in Hungary is 23% wheat, 26% grain maize, 14% sunflower, 7% barley, 5% rape and 7% fodder crops inside the arable lands (Hungarian Central Statistical Office [40]).

#### *2.3. Landscape Metrics*

The HEB database was applied to calculate size- and shape-related landscape metric parameters. Patch-level landscape indices were calculated for each LULC patch of the HEB database with the V-LATE 2 extension of Arc GIS 10.3 software [41]. Patch level metrics, created for individual land cover patches, characterize the spatial character and context of patches. These patch metrics serve primarily as the computational basis for developing a landscape metric. During our landscape metrics analyses, we calculated the following patch-level landscape metrics, which represent size and shape characteristics of land-cover patches (Table 1). The mean patch size (MPS) has been widely applied in landscape ecology, since it is commonly agreed that the occurrence and abundance of different species and species richness strongly correlates with the mean patch size. The shape complexity of individual LULC types was quantified by using landscape metrics (MFRACT). We applied the Shannon Diversity Index (SDI) to determine the landscape heterogenety [25]. We calculated these landscape indices (MPS, MFRACT, SDI) inside the 600 m radius buffer zones.


**Table 1.** Descriptions and calculations of the applied landscape indices [26,42,43].

#### *2.4. Statistical Analyses*

To understand the relationship between LULC types and skylark abundance, first we had to identify those LULC categories which are selected (used as habitat) by skylark or are avoided. We applied a preliminary test to identify the group of correlated land-cover and landscape index variables using variance inflation factors (VIFs), and the explanatory variables were not linearly related. VIF values were between 0 and 1.9, which shows that the multicollinearity is low between the variables (LULC types and indices). The arableland category was ignored from statistical analyses (model) because in Hungary and other European countries, the agricultural land is the matrix (dominant LULC type) in the landscape, so the proportion of this category shows strong autocorrelations with other LULC types. We used generalized linear models (GLM) to determine the impact of land cover and landscape structure (composition) on skylark abundance. We applied negative binomial models (link = log) to account the overdispersion of skylark-abundance data (tested by overdispersiontest function of AER package in R). Models with all possible combinations of explanatory variables were generated, and we established Akaike's information criterion to rank them with the "dredge" function from the MuMin package in R [44]. We used model averaging for competitive models (delta AICc < 2) to include uncertainty arising from the high number of candidate models (Table A3) [45]. The significance of the variables was estimated by the LmerTest package [46]. We constructed two groups from the LULC categories of the HEB database based on GLM results, namely, preferred (significant positive relation) and nonpreferred (significant negative relation) land-cover types. We analyzed the relationship between the landscape metrics of the preferred (as habitats) and nonpreferred land-cover types, and the skylark abundance data with negative binominal GLM and model averaging. In the next step in our investigation, we analyzed the shape and size characteristics of those LULC types which showed significant positive relation with skylark abundance. These land-cover types were separated into arable lands and

grasslands because we wanted to analyze the effect of arable land and grassland landscape metrics on the skylark population. In this model, the arable land category has been used. The distribution of landscape metric variables was not normal, so logarithmic transformation was used to normalize the data. These variables were in different dimensions, so we created a range function in R that transformed the variable values into a number between 0 and 1:

$$\text{Range function} = \frac{\mathbf{x} - \min(\mathbf{x}, na.rm = T)}{\max(\mathbf{x}, na.rm = T) - \min(\mathbf{x}, na.rm = T)},$$

where *Range function* is a number that describes the given number between 0 and 1, *na.rm = T* means that NA values were removed, *min* is the minimal value of the list, and *max* is the maximal value of the list. On the basis of the output of the statistical model, we could describe the optimal landscape configurations for this species.

#### *2.5. Model Validation*

We calculated the predicted marginal effects (ggeffects package in R) of the preferred land-cover types and their landscape metrics on the skylark population [47]. To validate our model, we set up a training and a testing group (66.6% and 33.3% proportion, respectively) with random sampling (sample.split function from caTools 1.17 package) on the basis of our dataset in R statistics software. We used the predict function from the car package to calculate the estimated skylark-abundance data. Model accuracy was measured by three indices: Spearman's rank correlation to show the relationship between observed and predicted values, mean absolute error to show the distance of the predicted values from the observed values [48], and mean absolute percentage error to show the percentage of error between observed and predicted values [49].

#### *2.6. Prediction of Skylark Population in Natura 2000 SPAs*

We could estimate skylark population density using the 600 m buffer areas and the HEB dataset. The centers of the buffer zones were in a regular grid (1200 × 1200 m) inside the Natura 2000 SPA dataset. We used the Natura 2000 SPA areas as the basis of our prediction site, because of the Eurasian skylark is a very common indicator species of agrarian landscapes (Natura 2000 Annex I. list). In Hungary the Natura 2000 SPA areas are typical agrarian landscapes which contain Urban areas (1.5%), Croplands (31.7%), Grasslands and other herbaceous vegetation (21.7%), Forest and woodlands (27.8%) and wetland and water surfaces (17.2%). Mowing of the grasslands inside the Natura 2000 sites is regulated by the law. The mowing machine should cut the grass 10 cm above the soil surface. Mowing should not begin before 1 of July, to protect the ground nesting birds. The number of the animals and the method (it is different based on the grassland type) are also regulated by the law. Prediction was performed based on the model results that analyzed the connection between the preferred area and the landscape metrics. Landscape indices were calculated inside the Natura 2000 SPAs. The estimated skylark population was calculated by the predict function in R software. Figure 3 shows the spatial distribution of the 600 m buffer zones.

**Figure 3.** Example the spatial distribution of the 600-m buffer zones inside a Natura 200 Spa protected area of Hungary.

#### **3. Results**

#### *3.1. Relationship between Land-Cover Proportions and Skylark Abundance*

Based on GLM results, we identified two main groups (classes) of the LULC categories of the HEB database. Preferred LULC categories that were considered the habitats of the Eurasian skylark because they showed significant positive relation with skylark abundance were those such as salt steppes and meadows, and closed grasslands in hills and mountains. The closed-grasslands LULC category showed the highest significant relation, thereby having the most important effect on skylark abundance. The arable-land LULC category is also a preferred category according to the literature [11,18,29,50]. The nonpreferred group (class) of LULC categories contains land-cover types with significant negative relations with skylark abundance: built-up land, green urban areas, complex cultivation patterns, forests, and wetlands and water surfaces. The complex-cultivation-pattern LULC category had the strongest negative association with the skylark population, followed by wetland and water surfaces, and green urban areas. The relative importance of the significant variables was 100% in all cases (Table 2).

models showing relative importance of each explanatory variable on Skylark abundance, estimated parameter values ± Standard deviation. (For detailed descriptions of the LULC categories see Table A1). **Relative**

**Table 2.** Summary table for LULC categories, which shows the GLM results after multimodel averaging of best candidate


Number of MMM observations (data pairs): 3049, \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001, Positive significant relation with skylark

abundance, Negative significant relation with skylark abundance, No significant relation with skylark abundance.

#### *3.2. Relationship between Landscape Structure (Compositon) and Skylark Abundance*

The landscape metrics that describe the shape and size characteristics of the preferred and nonpreferred LULC classes showed different directions of significant relation with skylark abundance (Table 3). The metrics that describe the shape complexity and size of the LULC patches of preferred LULC categories of the HEB database showed significant positive relations with skylark abundance. The shape complexity (MFRACT index) of the preferred LULC patches has stronger influence on the skylark abundance than the mean patch size (MPS). The shape complexity and size of the nonpreferred LULC categories had significant negative relation with skylark abundance in this case, MPS had higher association with skylark abundance. (Table 3). Land-cover heterogeneity, described with SDI, had a significant negative effect on skylark abundance, which showed that this species prefers a homogeneous landscape.

**Table 3.** Summary table for landscape metrics, which shows the GLM results after multimodel averaging of best candidate models showing relative importance of each explanatory variable on Skylark abundance, estimated parameter values ± Standard deviation.


Number of MMM observations (data pairs): 3049, \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001.

#### *3.3. Impact of Preferred Land-Cover Categories and Their Landscape Metrics*

Total grassland proportion had the highest association with skylark abundance, as shown in Table 2; the average size of arable-land patches (MPS) was more important from an abundance point view of this species than the mean patch size (MPS) of grassland patches. The complexity of grassland patches (MFRACT) had a significant positive association with skylark abundance, while the shape characteristics of arable land had no significant relationship with skylark abundance. The predicted marginal-effect graphs visualize the above-described connections between proportions of LULC categories, sizeand shape-related landscape indices, and the estimated population density changes of the

skylarks (Figure 4). According to the modeled population density changes, in the case of 100% grassland coverage of a hypothetical landscape, we could find about 4–6 skylark individuals/km\*. While the connection between the change in proportions of different land-cover types showed a near exponential curve, landscape metrics showed almost flat linear connections with estimated skylark abundance.

**Figure 4.** Predicted marginal effects between the skylark individuals / km2 proportions and landscape metrics of arable and grasslands. The confidence intervals (95%) of the prediction are shown between the dotted lines. ((**A**), Connection between the proportion of arable land and estimated population density of skylark, (**B**), Connection between the proportion of grassland and estimated population density of skylark, (**C**), Connection between the MPS of arable land and estimated population density of skylark, (**D**), Connection between the MPS of grassland and estimated population density of skylark, (**E**) Connection between the MFRACT of grassland and estimated population density of skylark).

On the basis of our results (Table A2), we could create an equation that describes and estimates the skylark population in a given landscape:

$$Skylark\_{\text{population}} = -3.24 + 1.29 \ast MPS\_{\text{variable land}} + 0.97 \ast MPS\_{\text{grossland}} + 0$$

$$0.63 \ast MFRACT\_{\text{grossland}} + 1.65 \ast Area\_{\text{arable land}} + 2.4 \ast Area\_{\text{grossland}}$$

where *Skylarkpopulation* is the skylark number density (individual/km\*), *MPSarable land* is the mean patch size of arable land, *MPSgrassland* is the mean patch size of grasslands, *MFRACTgrassland* is the mean fractal dimension of grasslands, *Areaarable land* is the proportion of arable land, and *Areagrassland* is the proportion of grasslands.

#### *3.4. Model Validation*

According the validation of our results, there was a significant Spearman's correlation between the observed and predicted skylark abundance values. Mean absolute error shows the distance between the predicted and observed abundance values of this species, which is +– 2.12. Mean absolute percentage error (MAPE) shows the prediction accuracy of the model in percentage; in this case, it was 37.77%. The accuracy of this model based on the MAPE was 62.23% (Table 4). If the model contains just the land cover types, the MEA is 2.95; MAPE is 46.56% and the Spearman correlation coefficient is 0.493.

**Table 4.** Summary table of the correlation and error indices, which show the accuracy of the predicted values, based on land cover types and land cover types + landscape indices.


#### *3.5. Prediction of Skylark Population of Natura 2000 Special Protection Areas of Hungary*

The spatial distribution of the predicted skylark population in each 600-m zone of the Natura 2000 SPAs of Hungary was very diverse (Figure 5). The total investigated Natura 2000 SPA was 13,514 km\*\* which cover the most valuable agroecosystems and rural landscapes of Hungary. Based on model prediction (predict function in R) inside these protected areas, approximately 23,746 skylark individuals were predicted. The density of this species is the highest in the agricultural-landscape-dominated areas of the great Hungarian plain.

**Figure 5.** Predicted Eurasian skylark population (individuals/km\*) in the 600 m buffer zones inside the Natura 2000 SPA area.

#### **4. Discussion**

There are several publications analyzing the relationship between skylark and LULC [4,51–54] in local small study areas, but our very detailed LULC dataset (HEB) offers a unique opportunity to obtain regional (country)-scale information about this relationship In our study, we considered both datasets describing proportions of LULC categories and landscape indices that describe the shape and size characteristics of preferred (habitat) and nonpreferred LULC categories. Based on our research findings, population density (individuals/km\*) could be estimated because there was a significant statistical relationship between proportions, the shape and size characteristics of different LULC types, and the abundance of this farmland bird. One new finding from our research is that, for the estimation of skylark population density, it is necessary to consider landscape indices together with the proportions of different LULC categories because shape (mean fractal dimension) and size (mean shape size) characteristics of these LULC categories also have significant association with skylark abundance. Based on our finding we have predicted the number skylarks inside the Natura 2000 SPA areas in Hungary.

#### *4.1. Impact of Proportions of LULC Categories on Skylark Abundance*

We could select two LULC groups (classes) from the land-cover types of a very detailed (20 × 20 m resolution) LULC map. Nonpreferred types had negative significant relation with skylark abundance. These were built-up and green urban areas, which negatively affected the population because of the lack of openness and the high proportion of constructed surfaces. Our findings are confirmed by other international publications [4,6,10]. The complex cultivation pattern land-cover type has negative significant relation with skylark data. Other authors underlines that the skylarks do not prefer heterogeneous agricultural lands because this rural landscape contains many different LULC patches, including also those that are not preferable to the skylark, like vineyards, fruit and berry plantations (because of its height, they obscure the view) [10,11,16,17]. Small parcels of, annual crops, city gardens pastures, fallow lands and/or permanent crops somewhere with scattered houses. Forest and wetland LULC categories are well-known nonpreferred land-cover types of the skylark. The skylark is a typical farmland bird; therefore, it is not a surprise that wetland areas, water bodies, and water courses are not suitable habitat types for this species. The main reason of the negative significant relation of the forest is the lack of openness, which is very important for the skylark [10,11,16,55]. In our research we were not take difference between the type of forests, because according to previous studies every types of forest areas are not habitats of this species.

In the estimation of skylark population density, the preferred land-cover types had higher weights (were more important) than those of the nonpreferred LULC categories. Arable land is a well-known habitat type of this farmland bird species according to the international literature [11,18,19,56,57]. Unfortunately, in Hungary is no available detailed country scale spatial statistical data about the cultivated crop types inside the arable lands (cropland) areas. According to the available most detailed Hungarian LULC dataset, the HEB dataset the 57% of Hungary is covered by agricultural fields and its 81% is arable land (Cropland). Grassland and pasture areas are also preferred LULC categories for skylark, [7,8,10,11,58–61]. The HEB dataset allow us to analyze the impact of different types of grassland on skylark abundance. We did not find significant statistical relations with open sand steppes and open rocky grasslands because the number of 600 m circle radius observation points of LULC categories have been low, and these landscape conditions (too-fragmented grassland areas with very short and very sparse vegetation) are not suitable for breeding skylarks [57,62]. There was a significant positive relation between skylark abundance, and the LULC categories of salt steppes and meadows, and closed grasslands. Each LULC category is suitable for nesting breeding skylarks because of the medium vegetation height and optimal proportion inside the 600 m radius circles. Our results are similar with those of others, who described strong relationship between closed grasslands and meadows and skylark abundance, the reason of this relation could be

the larger amount of food [63–66]. According to our findings for the prevention of the farmland bird habitats, the EU agri-environmental policy should pay more attention to the management of salt steppes and meadows, and closed grasslands. To increase the population density of skylark, the mean patch size and the proportion of these land cover types (compare to all) in the landscape should increase. In case of the protected grassland areas, one of the biggest ecological problem is the spontaneous spreading of the bush vegetation, which can reduce the skylark habitats. If we want to stop this process, and keep the openness of the landscapes, we should reduce the size and the shape of the bush and forest patches inside these grassland areas. Therefore, we must eradicate the spontaneously spread bush vegetation (which often full of invasive species) by the proper way of grazing or haymaking, the grasslands can keep its size, shape, and openness characteristics in the protected landscapes. This kind of management of protected areas can preserve not only the vegetation diversity of grasslands but it has also important key factor in the skylark habitat protection.

#### *4.2. Impact of Land-Cover Categories and Their Landscape Metrics*

The landscape metrics of the preferred LULC classes showed positive significant relation with skylark abundance, meaning that, if arable-land and grassland proportion and shape complexity was higher, then the skylark population would also be higher. The landscape metrics of the nonpreferred LULC classes showed negative significant relation with the skylark population, meaning that, in landscapes with small size and in compactshape nonpreferred LULC categories, skylark population density (abundance) would be higher.

LULC landscape heterogeneity has a negative effect on the skylark in this scale, where one land cover patch can contain more parcels. If landscape heterogeneity increases, the skylark population declines. This species prefer the homogenous LULC structures, which is in accordance with the results of other authors [10,11,16,17,62].

The grassland proportion had the highest association with the skylark population. This species usually nests and feeds in grasslands. The proportion of arable land has a high association with skylark abundance, but the level of its significance is lower. In the case of the MPS, the opposite phenomenon was observed: the MPS of arable lands (arable land patches of HEB) had a higher effect on skylark abundance than that of grassland. The skylark does not prefer small size arable lands (parcels) and grassland fields in that scale, where one arable land patch can contain more parcels [7,51,60,64]. According to Uuemaa et al. 2009 most bird species react more strongly to the composition land cover than to the configuration of landscapes [25]. Our results also show that the LULC proportions and mean patch sizes have stronger impacts on the abundance of this species, than the shape (fractal dimension index) characteristics of the habitat patches. The mean-absolutepercentage-error value (37.77%) was acceptable since, for a more precise prediction, we would have to use more variables (e.g., species and quantity of insects, used pesticides, parcel management) that are not accessible in country-scale analysis. We can determine that the landscape indices improved the model accuracy, based on the Table 4.

#### *4.3. Predicted Population Inside the Natura 2000 SPAs*

In Hungary, the latest estimated country-wide Eurasian skylark population is from 1999–2002. There is no spatially detailed population estimate. This study is the first estimate for Natura 2000 protected areas in Hungary. There some early 2000s studies about the skylark densities in Europe (Table 5).

The studies listed above do not use the shape and size related landscape indices for estimation of the skylark abundance (density). With the combination of the detailed point-based bird census data, detailed country-wide LULC dataset and landscape indices we can get a more precise prediction of skylark population. Our results are comparable with these previous estimations and the density values are similar [63,67–69].


**Table 5.** Summary table of studies, which predicted the Eurasian skylark density inside European study areas.

#### **5. Conclusions**

Landscape composition (proportions, and shape and size characteristics of LULC categories) has significant association with the skylark population. The salt steppes and meadows, and closed grassland serve as habitat for the Eurasian skylark. This study provides new information about the relationship between landscape metrics of the habitat types (shape and size characteristics of patches) and skylark abundance. Fractal dimension index, which describes the shape complexity of grassland patches has a positive impact on the skylark abundance, while the shape complexity of non-habitat types shows opposite relationships with the skylark density. We analyzed them together and could estimate the association of these landscape composition variables (proportions, shape and size characteristics of LULC classes) with skylark abundance. We could estimate skylark population density inside Natura 2000 SPAs in Hungary.

The outcomes of this study can be used for further land use planning, and the habitat design of Natura 2000 SPAs and other protected areas of the rural landscapes. According to our findings, inside the protected areas should increase the proportion, the average size and shape complexity of those LULC types (arable land, salt steppes and meadows, and closed grassland), which shows positive relations with the abundance data of skylark. It is feasible by stopping the spontaneous reforestation and eradicating the spontaneously spread vegetation (especially invasive bush species). The grazing or mowing, the protected grasslands can preserve the size, shape and openness characteristics of these skylark habitats. This kind of environmental management forms help to conserve the habitat types of skylark. The skylark is an area sensitive species and it is an indicator species of farmlands, so the shown methodology is adaptable for analyzing the impact of landscape composition on other farmland bird populations [56,70–72]. The skylark is considered as indicator for monitoring of agricultural landscapes, because its abundance shows strong relationships with other farmland bird species [73].

**Author Contributions:** Conceptualization, N.C. and P.S.; methodology, N.C.; validation, N.C.; writing—original draft preparation, N.C. and P.S.; writing—review and editing, N.C. and P.S.; visualization, N.C.; supervision, P.S.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the "UNKP-20-3-SZTE-515 New national excellence program of the Ministry for Innovation and Technology".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We thank the anonymous referees for their valuable recommendations and suggestions. We also thank Károly Nagy and MME (Bird Monitoring Centre, MME/BirdLife Hungary) for data provision.

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

#### **Appendix A**


**Table A1.** The LULC categories of the Hungarian Ecosystem Basemap, and the investigated LULC categories.

**Table A2.** Summary table for landscape metrics and LULC categories, which shows the GLM results after multimodel averaging of best candidate models showing relative importance of each explanatory variable on Skylark abundance, estimated parameter values ± Standard deviation. MPS is Mean Patch Size and MFRACT is Mean Fractal dimension.


Number of MMM observations (data pairs): 1897, \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001.

**Table A3.** Summary table of component models from model averaging.


1 Built-up, 2 Green urban areas, 3 Permanent crops, 4 Complex cultivation pattern, 5 Open sand steppes, 6 Salt steppes and meadows, 7 Open rocky grasslands, 8 Closed grasslands in hills and mountains or on cohesive soil, 9 Other herbaceous vegetation, 10 Forest, 11 Wetlands and water surfaces.

#### **References**


### *Article* **Structural Variations in the Composition of Land Funds at Regional Scales across Russia**

#### **Vasilii Erokhin 1, Tianming Gao 1,\* and Anna Ivolga <sup>2</sup>**


Received: 12 May 2020; Accepted: 17 June 2020; Published: 17 June 2020

**Abstract:** In recent decades, Russia has experienced substantial transformations in agricultural land tenure. Post-Soviet reforms have shaped land distribution patterns but the impacts of these on agricultural use of land remain under-investigated. On a regional scale, there is still a knowledge gap in terms of knowing to what extent the variations in the compositions of agricultural land funds may be explained by changes in the acreage of other land categories. Using a case analysis of 82 of Russia's territories from 2010 to 2018, the authors attempted to study the structural variations by picturing the compositions of regional land funds and mapping agricultural land distributions based on ranking "land activity". Correlation analysis of centered log-ratio transformed compositional data revealed that in agriculture-oriented regions, the proportion of cropland was depressed by agriculture-to-urban and agriculture-to-industry land loss. In urbanized territories, the compositions of agricultural land funds were predominantly affected by changes in the acreage of industrial, transportation, and communication lands. In underpopulated territories in the north and far east of Russia, the acreages of cropland and perennial planting were strongly correlated with those of disturbed and barren lands. As the first attempt at such analysis in Russia, the conversion of cadastral classification data into land-rating values enabled the identification of region-to-region mismatches between the cadaster-based mapping and ranking-based distribution of agricultural lands.

**Keywords:** agricultural land; cropland; land category; land fund; territory; Russia

#### **1. Introduction**

Structural alterations in land use have been intrinsically associated with a growing demand for food [1,2]. Increasingly, contemporary processes of progressing urbanization and industrialization have been aggravating the conflicts between different functional land types [3]. As land systems represent a critical intersection between economic and ecological systems [4], land distribution patterns are becoming more vulnerable to a variety of environmental and social issues. Out of the one-fifth of the world's total land surface, which is potentially suitable for crop production, more than half is already actively cultivated [5]. Further agricultural expansion is hampered by natural and geographical factors [6], pervasive land-use change impacts [7], high economic costs [8], and infrastructure constraints [9]. At the same time, according to DeFries et al. [10], Ajani [11], and Lambin [12], agricultural production tends to face increasing competition for land with other types of land use. Over recent decades, many scholars and practitioners, including Platt [13], Briggs and Yurman [14], Vining et al. [15], and Sioen et al. [16], among others, have been reporting the irreversible removal of substantial areas of land previously used for agriculture to urban, industrial, infrastructure, and other types of use instead. Urbanization and industrialization intensify competition between agricultural and non-agricultural land-use practices [17]. Along with industrial development and urban sprawl, there are significant alterations of land use far beyond city limits that result in arable land loss [18].

Generally, at a regional scale, agricultural lands do not strictly compete with other categories for the same land areas due to the specific climate, soil, and topographical requirements for farming. However, in land-abundant and climate-diverse countries, the geographical distribution of agricultural land use tends to adjust to better match land quality [19]. Russia is aa good example of aa country that can be used to demonstrate this fact. Agriculture abandonment in vast northern and eastern areas has occurred in parallel with a concentration of intensive agriculture in fertile lands in the southern, western, and central regions of the country. In Russia, agricultural lands only represent 12.96% of the total national land fund (cropland at 7.16%, rangeland at 3.99%, hayfields at 1.40%, fallow at 0.28%, and perennial plantings at 0.11%). Per-territory concentrations of agricultural land vary from 75.32% in the Southern Federal District and 70.96% in the North Caucasian Federal District to only 4.05% in the Northwestern Federal District and 1.30% in the Far Eastern Federal District.

We clarified the definitions of the main terms used in this study as follows:


The disproportions of agricultural land distribution are, to some extent, caused by economic factors, not only geographic and natural conditions. Similar to most post-socialist countries, Russia has experienced dramatic changes in land ownership and land tenure since the early 1990s. Among the principal transformations, Lerman and Shagaida [20] have outlined the privatization of agricultural land, rights to agricultural land for individual landowners, and the removal of prohibitions on buying and selling land. The land market has responded positively to the liberalization with an increase in transactions between individual landowners [20]. However, the domination of shared and joint land ownership has weakened the role of the state in controlling land use [21] and has increased the fragmentation of public land property into many scattered units [22]. Almost twelve million land shares (certificates) were distributed between rural individuals and former employees of collective and state farms [23]. According to Trukhachev et al. [23], Lerman and Shagaida [20], Rozhkov [24], and Visser et al. [25], land reform in Russia has significantly contributed to structural variations in the composition of land funds. The proportion of agricultural land in the total land fund has been declining due to a loss of arable land, particularly in the vast areas of the Far Eastern Federal District and the Siberian Federal District [26]. From 1990 to 2000, the rate of land abandonment in Russia was above 30%, one of the highest among the economies in transition [27]. Milanova [28] reported a decrease in the cropped area for all crops during the 1990s due to the changes in land tenure and stagnation of the agricultural sector. A drastic decline in livestock production resulted in a reduction of hayfields and rangelands. Vast areas of arable land were abandoned due to land degradation. In some territories in the central, northern, and eastern parts of the country, humus content dropped by 50%. Prishchepov et al. [29] revealed the correlation between the spatial distribution of abandoned croplands and natural factors, such as inadequate precipitation and shorter growing periods, in both Siberia and eastern parts of the country. As many farms were situated in the boreal zone, some of the abandoned lands have experienced shrub and tree encroachment [30].

Many experts report an aggravated environmental degradation of agricultural lands due to over-exploitation [31,32]. The changes in land cover and land use in forest-steppe and steppe vegetation zones (agriculture-oriented territories of southern Russia, the European center, and southern parts of Ural and Siberia) have been driven by extensive farming. Milanova [28] and Milanova et al. [33] reported that up to 90% of lands in some territories were converted to crop production. However, where environmental concerns of land use are mentioned in either federal or regional legislation, they predominantly relate to reducing industrial emissions or waste disposal in urban and suburban areas, not to agricultural land use [34]. Over 40 million hectares of cropland is now abandoned in Russia, and another 58 million is eroded. Land degradation, along with desertification due to irrational land use, poses serious environmental, economic, and social threats in the long-term. Griewald et al. [34] argued that the land use context in Russia did not support a transition towards sustainable land management, i.e., a "use of land resources, including soils, water, animals, and plants, for the production of goods to meet changing human needs, while simultaneously ensuring the long-term productive potential of these resources and the maintenance of their environmental functions" [35]. The urban expansion causes shrinkage of arable and other categories of agricultural land [36], which are transferred to various non-agricultural types of land use. A considerable amount of agricultural land loss due to urbanization and industrialization takes place on fertile soil [37] and irrigated lands [38]. In return, the increase in agricultural land acreage occurs on soils that are lower in terms of their fertility. Prishchepov et al. [39], Brueckner [40], and Brown et al. [41] raised concern over the growing concentration of arable land in smaller and more fragmented locations in proximity to urban and industrialized areas. Erma et al. [42] reported many cases where residential settlements occupied agricultural land in southern and central parts of the country, which are known as the breadbasket regions of Russia.

With increased variability in the composition of land funds, a reliance on research in this area has become more critical. In a series of empirical studies, many authors, including Verburg et al. [43], Van Doorn and Bakker [44], Nainggolan et al. [45], and Diogo and Koomen [46], among others, have attempted to construct hypotheses about the relationship between proximate driving forces and agricultural land-use patterns. The problem is that the established hypotheses do not adequately explain the causality between land-use processes and the compositions of land funds at different regional scales. In transition economies, including Russia, where land reforms have dramatically changed the distribution of the land inventory in recent decades [42], variations in agricultural lands due to the pressure of non-agricultural land use have remained under-investigated. The composition of agricultural land funds has commonly been considered out of a non-agricultural context [47,48], instead of exploring the interactions between the proportions of agricultural, urban, infrastructure, and industrial lands. Most of the studies have applied a proportion of agricultural land in a land fund as a core territorial specification without further testing for alternative non-agricultural land use variables [4]. Therefore, in regional studies, a knowledge gap has emerged in terms of how the variations in the compositions of land funds may be tracked with an aim to optimize agricultural land use. A more explicit focus on the relationships between land categories is required to be able to explain and predict land system dynamics in diverse locations [49]. With this background, in the case of Russia, this study aimed to contribute to the body of knowledge on regional scale land uses by identifying structural variations in the compositions of territory land funds and revealing the interdependencies between the proportions of agricultural, on the one side, and urban, industrial, and other types of land on the other.

#### **2. Materials and Methods**

This study was a quantitative study that was performed based on the data obtained from land registers from 82 out of 85 of the administrative entities of Russia (further detailed in Section 2.6). Russian public statistics report thirteen land categories within land funds, including five agricultural and eight non-agricultural ones. As we aimed to study structural variations in the compositions of land funds by identifying the changes in the proportions of different lands, all thirteen land categories were considered here (the definitions are given in Section 2.1). The overarching methods adopted in this study included a ranking of the territories on the degree of agricultural land activity (see Sections 2.2 and 2.3), centered log-ratio transformation of compositional land share data to an unconstrained space

(Section 2.4), and correlation analysis to reveal the variations in the proportion of land categories within the groups of territories (Section 2.5). In total, the study algorithm followed five stages (Table 1), which are further addressed in Sections 2.1–2.5 of the paper.


Source: Authors' development.

#### *2.1. Stage 1: Land Categories*

As the structural features of land classification frameworks largely depend on the purpose of classification [50], various country specific approaches exist to categorize agriculture and other types of land. In Russia, Shagaida [51], Nosov [52], and Macht et al. [53] have contributed to the identification of various categories of agricultural lands. The majority of the studies, however, have paid inadequate attention to revealing variations in land fund compositions due to the specific needs for farming, residential construction, or industrial and infrastructure development in particular locations. For instance, Zhang et al. [3] applied an ecological-living-production classification system, to demonstrate the distribution of agricultural land across arable land, pastures, timberland, aquaculture land, and orchards, but they did not reveal the variations in the spatial concentration of particular land categories. Loshakov [54] developed an approach for the categorization of agricultural lands based on the productive qualities of soils but did not consider mismatches between agro-climatic zoning and land registers.

While the adjustments to land classification systems may be useful in achieving some specific technical, geographical, environmental, or economic goals, there are situations in which various existing approaches should be merged [55]. Many of the systems have a limitation in their ability to demonstrate the interrelationships between the categories of land cadasters for agricultural production. In general, classification concepts do not correctly emphasize per-category changes in the composition of a land fund. This is also one of the inherent vices of state statistics reporting on land fund structures in many countries. Notably, the Federal State Statistics Service of the Russian Federation (Rosstat) generalizes land into three broad categories (namely, agricultural, woodlands, and water reserve lands) [56]. Separate forms also report urban lands and lands for industrial, transportation, and communication infrastructure purposes; however, these forms exist at a national scale, not a regional scale. More detailed classification for five categories of agricultural land (croplands, hayfields, rangelands, perennial plantings, and fallows) is available in agricultural census report forms [57]. However, since the agricultural census is conducted decennially, intercategory variations cannot be effectively tracked on an annual basis.

One of the possible solutions to this discontinuity problem is to supplement census data with operative land cadaster information [58]. In Russia, the Federal Service for State Registration, Cadastre, and Cartography (Rosreestr) continually monitors land fund compositions per territories across seven categories of land, including agricultural land, residential land, industrial land, specially protected territories, woodlands, water fund lands, and reserve lands [59]. Among several classification schemes used by Rosreestr, one breaks agricultural lands into five categories, similar to Rosstat's decennial census, but instead on an annual basis. The usage of this data may allow the creation of a better time-sensitive model to represent changes in the proportion of land categories within different regions.

In this study, simple classifications determining the allocation of land between agriculture, urban, and nature were merged with more comprehensive ones, in which cadaster synergies could be detailed for a wider range of agricultural, industrial, urban and built-up, forest, and water reserve lands. The array included the categories of urban and infrastructure lands (obtained from separate sections of Rosstat's reports), as well as wetlands, disturbed lands, and barren lands (all reported by Rosreestr's alternative classification of utilized lands). In total, the authors' model merged thirteen land categories, including five agricultural (*L*(1–5)) and eight non-agricultural (*L*(6–13)) categories (Table 2). As reported by Rosstat [56,57] and Rosreestr [59], the categories were mutually exclusive and exhaustive. That is, each location within the *Tj* territory could be classified into one and only one *Li* category.



Source: Authors' development based on Rosstat [56,57] and Rosreestr [59].

#### *2.2. Stage 2: Composition of Land Funds*

As the keynote idea is to reveal the variations in the compositions of the land funds across diverse territories, a kind of assessment scale should be applied. There have been many attempts to find a reliable approach for the conversion of cadastral classification data into land-rating values. Land classification systems based on rankings have been in use since the 1980s when Wright et al. [60] and Cocks et al. [61] first applied simple additive linear models of factor weights to the evaluation of land utility for crop production. In the realm of building a relevant ranking framework, one of the major challenges is determining how to align categorization (public statistics) with functional scales. In agriculture, variations between the proportions of lands are hard to identify [62] and thus cannot be effectively linked with territory fragmentations of agricultural production [63]. The immediacy of the problem was convincingly demonstrated by Grˇcman et al. [64], who found the difference between land-rating values based on precise calculations and those based on official information (specifically, for agricultural land with lower production potential).

Another challenge is that the ranking systems are not comparable and, therefore, inapplicable across a variety of agricultural and non-agricultural lands [65,66], and even across croplands, fallows, and pastures [67,68]. There have been attempts to overcome this problem by finding an integral parameter that would allow the adjustment of agricultural- and non-agricultural-oriented ranking systems to be comparable. In terms of land fund compositions, one of the most promising foundations of ranking is the contribution of a land category to the total land acreage per territory [69] (Equation (1)). The applicability of this parameter for building category-based land assessment frameworks was successfully tested by Mazurkin and Mihailova [70], Buckett [71], Artamonova et al. [72], Stupen et al. [73], Shishkina et al. [74], and Yerseitova et al. [75].

$$A\_{j\text{Li}} = \frac{S\_{j\text{Li}}}{S\_j} \tag{1}$$

where *AjLi* = share of land category *Li* in the land fund in territory *Tj*; *SjLi* = area of *Li* in territory *Tj*; *Sj* = total land acreage of territory *Tj*.

The shares of the *L*(1–13) land categories in the land funds were computed across *T*(1–82) territories (Appendix B, Tables A9–A16).

#### *2.3. Stage 3: Agricultural Land Activity*

Further, the *AjLi* values are ranked across the arrays of *Li* land categories and *Tj* territories to calculate a parameter of land activity. Agricultural land activity is a degree of orientation of a land fund composition toward an agricultural type of land use. It is an indicator of how a proportion of *L*(1–5) to *L*(6–13) serves the purpose of agricultural production in particular geographic and economic conditions at a regional scale. Land activity is a score of a *Tj* territory, obtained based on the proportions of various land categories within a land fund. Higher contributions of *L*(1–5) to total acreage result in higher agricultural land activity scores. The activity-rank correspondence is straightforward, where the higher is *AjLi* value, the higher is *Rji* score. A high rank demonstrates an orientation of land fund composition towards agricultural specialization. Since the prevalence of non-agricultural lands is considered as a spatial constraint for the allocation of agricultural land uses, higher proportions of *L*(6–13) within a land fund result in lower agricultural land activity scores. For these land categories, the activity-rank relationship is inverse, where the higher is *AjLi* value, the lower is the *Rji* score. For *j* territories included in the study, the *R* interval was [0; *j* − 1]. In our model, as *AjL*(1–5) tended to 1, *R* tended to (*j* − 1), while as *AjL*(6–13) tended to 1, *R* tended to 0.

Then, we assessed the significance of derived estimates. *Rji* scores were used to identify the quartiles of *AjLi* (Figure 1). The [- *Rjmin*; - *Rjmax*] interval was divided into the quartiles by finding the *n* multiplier, where *n* = - *Rjmax*− - *Rjmin* <sup>4</sup> .

**Figure 1.** Scale to classify *Tj* territories on the degree of agricultural land activity. Source: Authors' development.

The quartile-based approach was used by Mazurkin [69] for the ranking of territories based on absolute values of land activity parameters. It also agrees with Kotykova et al. [76] and Zhildikbaeva et al. [77], who compared the deviations of land category estimates from their highest level on a territory-by-territory basis. In this study, such a method for the classification of rankings allowed consideration of the information in the percentage areas measured for each *Li* in each *Tj*.

#### *2.4. Stage 4: Revealing Structural Variations of Agricultural and Non-Agricultural Land Categories*

Since the early years of Russia's land reform, structural variations in the compositions of land funds have progressed in response to socioeconomic and anthropogenic processes. To identify these variations between various land categories across four types of territories, this study employed factor analysis. It enables the transformation of land fund data into meaningful information [43,78] and revelation of variations in the structure of the use of territory land funds. According to Alcamo et al. [79] and Lavalle et al. [80], the integration of proximate and underlying factors may capture both the spatial distribution and the variety of land categories claimed for different land-based activities. The employment of factor analysis tools at a regional scale by Bakker et al. [81], Van Doorn and Bakker [44], and Hatna and Bakker [82] demonstrates the appropriateness of the method for cross-territory comparisons.

Among numerous factor analysis approaches, correlation analysis is one of the most suitable approaches to reveal variations in land fund compositions [83,84]. Since the *AjLi* data are compositional, i.e., they add up to a constant value of 1 or 100% of a land fund, they need a special treatment prior to correlation analysis [85]. Aitchison [86] named land fund compositions among the typical datasets associated with challenging problems in compositional data analysis. In a compositional vector that consists of several parts summing up to a constant, the relevant information is contained only in the ratios between these parts [87] (Equation (2)).

$$\mathbf{x} = (\mathbf{x}\_1, \dots, \mathbf{x}\_D)^t, \mathbf{x}\_i > 0,\\ \dot{\mathbf{z}} = 1, \dots, D,\\ \sum\_{i=1}^D \mathbf{x}\_i = k \tag{2}$$

where *D* = number of compositions, and *k* = a positive constant value, i.e., the sum of *D* compositions.

If correlation analysis is applied directly to the *AjLi* data, this can give misleading results [88] and form undesirable properties, like scale dependence [89]. The best way to analyze data with constant sum constraints is by first transforming them into an unconstrained space [88], where standard data analysis tools can then be employed [90]. Several log-ratio transformations have been introduced by Aitchison [89,91], Pawlowsky-Glahn et al. [87,90], Filzmoser and Hron [85], Long and Wang [92], and Van den Boogaart and Tolosana-Delgado [93]. Commonly used methods include using the additive log-ratio (alr), isometric log-ratio (ilr), and centered log-ratio (clr). Additive log-ratio transformation is based on log-ratios to a single reference variable. It is the simplest way to transform compositional data. However, it does not preserve distances between variables; i.e., it is not isometric [85]. Isometric log-ratio transformation is built on the choice of an orthonormal basis and thus solves the isometry problem. However, according to Egozcue et al. [94] and Egozcue and Pawlowsky-Glahn [95], base compositional parts are only related to isometric log-ratio transformed variables through non-linear

functions. In our case, this meant that the computed correlations between the proportions of land categories could not be interpreted in the sense of the *AjLi* data.

For this study, we employed centered log-ratio transformation (Equation (3)). Distinct from the additive log-ratio method, the centered log-ratio method is based on the geometric mean of all variables. It allows for the selection of a ratio variable to be avoided [85]. In contrast with the isometric log-ratio method, the centered log-ratio method simplifies the interpretation of the transformed variables because one could think of them in terms of the original variables [85,96].

$$y = [y\_1, \dots, y\_D] = \left| \ln \frac{\mathbf{x}\_1}{\sqrt[D]{\prod\_{i=1}^D \mathbf{x}\_i}}, \dots, \ln \frac{\mathbf{x}\_D}{\sqrt[D]{\prod\_{i=1}^D \mathbf{x}\_i}} \right| \tag{3}$$

where *x* = *AjLi* share of land category *Li* in the land fund in territory *Tj*; *y* = transformed *AjLi* compositions *ATRjLi*; *D* = number of compositions, i.e., *Li* land categories.

The *AjLi* compositions were transformed into *ATRjLi* data across all *Tj* territories using CoDaPack. This open-access software is one of the easiest-to-use applications that is commonly employed for compositional data transformation (for instance, see Thió-Henestrosa and Martín-Fernández [97], Egozcue and Pawlowsky-Glahn [98], and Muriithi [99]). The centered log-ratio-transformed data that were obtained were standard multivariate data that enabled us to use correlation analysis. Correlation matrices were built separately for the four groups of territories earlier ranked by the type of agricultural land activity. Correlation analysis was carried out here using the Excel Data Analysis ToolPak.

#### *2.5. Stage 5: Significance of Correlations*

When conducting correlation analysis for land systems, most scholars have faced a challenge similar to what we outlined earlier concerning ranking scales, namely, determining the significance of synergies between variables. Among various methods, the coefficient of correlation variance seems to be the most appropriate for dealing with interdependent multitudes of land categories [69,70] (Equation (4)).

$$C\_{cv} = \frac{\sum ATR\_{jLi}}{ATR\_{\text{max}} \times N\_L \times N\_T} \tag{4}$$

where *Ccv* = coefficient of correlation variance; - *ATRjLi* = sum of transformed *AjLi* values of *Li* land categories in *Tj* territories in the group; *ATRmax* = the highest value of *ATRjLi* in the group; *NL* = number of land categories in the array; *NT* = number of territories in the array.

The *Ccv* value was applied across four correlation matrices (types of land activity) to remove weak interdependencies and reveal strong synergies between the proportions of the *L*(1–5) and *L*(6–13) land categories in a land fund.

#### *2.6. Territories and Data*

Russia is a federation comprised of 85 administrative entities, or territories, as defined in the Section 1. Our study included 82 of them (mapped in Figures 2 and 3). The three municipal areas of Moscow, Saint Petersburg, and Sevastopol were excluded from the array as they are areas in which the proportion of agricultural land in the territory land fund is of negligible importance. For each territory, land cadaster data were derived from the annual reports from Rosreestr [59] and Rosstat [56,57] during 2010–2018. In Russia, these data are reported across thirteen land categories in thousand hectares. Appendix A summarizes the data of the total acreages of the territories included in the study, along with the acreages of the thirteen land categories. The study was built on the mean acreages of *L*(1–13) land categories during 2010–2018 (Appendix A, Tables A1–A8). The proportions of the *L*(1–13) land categories in regional land funds across *T*(1–82) territories are provided as percentages in Appendix B, Tables A9–A16. The variations in the proportions are provided as differences between 2010 and 2018 in Appendix B, Tables A9–A16. The consideration of the Republic of Crimea as a part of the array was

determined by the current position of the territory as being de-facto controlled by Russia. In no way, these results reflect the authors' attitude to the international status of the area. For the Republic of Crimea, we used the mean data of the land acreage and land categories' proportions from 2015–2018.

**Figure 2.** Spatial distribution of agricultural lands in Russia. Note: 1 = Belgorod; 2 = Bryansk; 3 = Vladimir; 4 = Voronezh; 5 = Ivanovo; 6 = Kaluga; 7 = Kostroma; 8 = Kursk; 9 = Lipetsk; 10 = Moscow Oblast; 11 = Orel; 12 = Ryazan; 13 = Smolensk; 14 = Tambov; 15 = Tver; 16 = Tula; 17 = Yaroslavl; 18 = Karelia; 19 = Komi; 20 = Arkhangelsk; 21 = Vologda; 22 = Kaliningrad; 23 = Leningrad; 24 = Murmansk; 25 = Novgorod; 26 = Pskov; 27 = Nenets; 28 = Adygeya; 29 = Kalmykia; 30 = Crimea; 31 = Krasnodar; 32 = Astrakhan; 33 = Volgograd; 34 = Rostov; 35 = Dagestan; 36 = Ingushetia; 37 = Kabardino-Balkaria; 38 = Karachaevo-Cherkessia; 39 = North Osetia-Alania; 40 = Chechnya; 41 = Stavropol; 42 = Bashkortostan; 43 = Mari El; 44 = Mordovia; 45 = Tatarstan; 46 = Udmurtia; 47 = Chuvashia; 48 = Perm; 49 = Kirov; 50 = Nizhny Novgorod; 51 = Orenburg; 52 = Penza; 53 = Samara; 54 = Saratov; 55 = Ulyanovsk; 56 = Kurgan; 57 = Sverdlovsk; 58 = Tyumen; 59 = Chelyabinsk; 60 = Khanty-Mansi; 61 = Yamal-Nenets; 62 = Altay Republic; 63 = Buryatia; 64 = Tyva; 65 = Khakasia; 66 = Altay; 67 = Zabaikalsk; 68 = Krasnoyarsk; 69 = Irkutsk; 70 = Kemerovo; 71 = Novosibirsk; 72 = Omsk; 73 = Tomsk; 74 = Sakha Yakutia; 75 = Kamchatka; 76 = Primorye; 77 = Khabarovsk; 78 = Amur; 79 = Magadan; 80 = Sakhalin; 81 = Jewish AO; 82 = Chukotka. The Republic of Crimea was included in the study due to its current position as a territory under the de-facto control of Russia. This in no way reflects the authors' attitude to the international status of the area. Source: Authors' development.

**Figure 3.** Russian territories: types of agricultural land activity. Note: 1 = Belgorod; 2 = Bryansk; 3 = Vladimir; 4 = Voronezh; 5 = Ivanovo; 6 = Kaluga; 7 = Kostroma; 8 = Kursk; 9 = Lipetsk; 10 = Moscow Oblast; 11 = Orel; 12 = Ryazan; 13 = Smolensk; 14 = Tambov; 15 = Tver; 16 = Tula; 17 = Yaroslavl; 18 = Karelia; 19 = Komi; 20 = Arkhangelsk; 21 = Vologda; 22 = Kaliningrad; 23 = Leningrad; 24 = Murmansk; 25 = Novgorod; 26 = Pskov; 27 = Nenets; 28 = Adygeya; 29 = Kalmykia; 30 = Crimea; 31 = Krasnodar; 32 = Astrakhan; 33 = Volgograd; 34 = Rostov; 35 = Dagestan; 36 = Ingushetia; 37 = Kabardino-Balkaria; 38 = Karachaevo-Cherkessia; 39 = North Osetia-Alania; 40 = Chechnya; 41 = Stavropol; 42 = Bashkortostan; 43 = Mari El; 44 = Mordovia; 45 = Tatarstan; 46 = Udmurtia; 47 = Chuvashia; 48 = Perm; 49 = Kirov; 50 = Nizhny Novgorod; 51 = Orenburg; 52 = Penza; 53 = Samara; 54 = Saratov; 55 = Ulyanovsk; 56 = Kurgan; 57 = Sverdlovsk; 58 = Tyumen; 59 = Chelyabinsk; 60 = Khanty-Mansi; 61 = Yamal-Nenets; 62 = Altay Republic; 63 = Buryatia; 64 = Tyva; 65 = Khakasia; 66 = Altay; 67 = Zabaikalsk; 68 = Krasnoyarsk; 69 = Irkutsk; 70 = Kemerovo; 71 = Novosibirsk; 72 = Omsk; 73 = Tomsk; 74 = Sakha Yakutia; 75 = Kamchatka; 76 = Primorye; 77 = Khabarovsk; 78 = Amur; 79 = Magadan; 80 = Sakhalin; 81 = Jewish AO; 82 = Chukotka. The Republic of Crimea was included in the study due to its current position as a territory under the de-facto control of Russia. This in no way reflects the authors' attitude to the international status of the area. Source: Authors' development.

#### **3. Results**

#### *3.1. Composition of Land Funds*

The analysis of land cadaster data across Russia's *Tj* territories (Appendix B, Tables A9–A16) allowed the discovery of a distinct regularity in the spatial distribution of agricultural lands. In southern and central parts of the country (green belt between 45◦ and 55◦ north latitude), croplands prevailed in the composition of the land funds (Figure 2). In the mountainous areas of North Caucasus, the blue belt comprised the territories where rangelands and other agricultural lands predominated. In most of the northern and eastern regions, the land funds were comprised of non-agricultural lands with a minor proportion of cropland.

#### *3.2. Agricultural Land Activity*

The ranking of Russia's territories on a parameter of agricultural land activity resulted in higher scores for the southern and central parts of the country than for Siberia and the Far East (Appendix C, Tables A17–A24). Concurrently, some less apparent findings were yielded (Appendix D, Table A25).

First, in the Southern Federal District, an agricultural granary for the country, the land fund composition was less agriculture-oriented compared to the Central and Volga districts and some territories of Siberia. Specifically, for Krasnodar and Rostov, two green belt territories with a considerable proportion of cropland in the structure of the land fund, the -*Rji* values were well below the district average. In some territories in the south and center, high ranks of cropland and rangeland were negated by low ranks for barren lands, water reserve lands, residential, industrial, and infrastructure lands.

Second, in the Siberian Federal District, the -*Rji* values nearly reached those values of the central and southern districts due to the high scores of hayfields in Omsk and Novosibirsk. The green belt by Altay was rated high for the proportion of cropland and other agricultural lands in the composition of the land fund.

Third, the yellow and red belts in the Far East feature the least agriculture-oriented macroregion in Russia. In Chukotka, Magadan, and Sakhalin, where woodlands and wetlands dominate the composition of the land fund, the agricultural land categories were ranked the lowest among the 82 territories examined here. However, in Primorye, Khabarovsk, Amur, and Jewish Autonomous Oblast, fallows, hayfields, and rangelands received high scores.

Following the obtained ranks, four *Rj* intervals were identified, each of which included *Tj* territories according to the degrees of agricultural land activity. The grouping reproduced the earlier revealed belt-like distribution of agricultural land, but with a modified configuration instead (Figure 3).

Generally, while the green belt shrank and shifted eastward, the blue one expanded and spread north of the 55◦ latitude mark. In some of the previously yellow belt territories of the Northwestern, Central, and Volga districts, perennial plantings and hayfields were ranked high enough to include those regions as type II regions. In Siberia, the green belt included Omsk and Novosibirsk due to the high rank of hayfields and the low rank of disturbed and barren lands. The blue belt stretched from Ural (Tyumen and Chelyabinsk) to Siberia (Tomsk, Khakasia, Tyva, and the Altay Republic) and farther to the Far East (Zabaikalsk and the Jewish Autonomous Oblast). In the south, the substantial activity of residential, industrial, transportation, communication, and disturbed lands downgraded Krasnodar to type III and Rostov and Crimea to type II. Ingushetia, Kabardino-Balkaria, Karachaevo-Cherkessia, and Dagestan, on the contrary, broke forth to the green belt due to high scores of perennial plantings and rangeland and low activity of wetlands, disturbed lands, and water reserve lands.

#### *3.3. Correlation Analysis*

In type I territories, the variations in the compositions of agricultural lands correlated with the changes in the acreage of non-agricultural land for infrastructure, primarily transportation and communication (the strongest correlation with cropland, perennial plantings, and hayfields) (Table 3). Strong correlations were also revealed between the proportions of croplands and fallows, on one side, and those of woodland and barren land on the other. The share of rangeland in the land fund was strongly correlated with that of barren land.


**Table 3.** Correlation matrix for type I territories.

Note: *ATRLi* = centered log-ratio-transformed data: *ATRL*<sup>1</sup> = cropland; *ATRL*<sup>2</sup> = fallow; *ATRL*<sup>3</sup> = perennial plantings; *ATRL*<sup>4</sup> = hayfields; *ATRL*<sup>5</sup> = rangeland; *ATRL*<sup>6</sup> = woodlands; *ATRL*<sup>7</sup> = forest range; *ATRL*<sup>8</sup> = water reserve lands; *ATRL*<sup>9</sup> = residential and industrial lands; *ATRL*<sup>10</sup> = lands under transportation and communication infrastructure; *ATRL*<sup>11</sup> = wetlands; *ATRL*<sup>12</sup> = disturbed lands; *ATRL*<sup>13</sup> = barren; bold denotes a strong correlation, *CATRli* > *Ccv* (0.7022 for type I territories). Source: Authors' development.

Similar to type I, in the type II group, a strong correlation was found between the shares of cropland and perennial plantings and those of lands for transportation and communication infrastructure (Table 4). Besides, since the blue belt predominantly was comprised of densely populated territories, there was a correlation between the shares of croplands and residential lands. In many type II territories, the contribution of woodlands and other forest ranges to the structure of the land fund was essential. This fact might explain the high correlation between the composition of agricultural lands and woodlands. In the south, where the climate and soil favor the development of horticulture and viniculture (i.e., in Crimea, Adygeya, and Rostov), *Ccv* emphasized a strong correlation between the proportions of perennial plantings and croplands within the agricultural land categories.

**Table 4.** Correlation matrix for type II territories.


Note: *ATRLi* = centered log-ratio-transformed data: *ATRL*<sup>1</sup> = cropland; *ATRL*<sup>2</sup> = fallow; *ATRL*<sup>3</sup> = perennial plantings; *ATRL*<sup>4</sup> = hayfields; *ATRL*<sup>5</sup> = rangeland; *ATRL*<sup>6</sup> = woodlands; *ATRL*<sup>7</sup> = forest range; *ATRL*<sup>8</sup> = water reserve lands; *ATRL*<sup>9</sup> = residential and industrial lands; *ATRL*<sup>10</sup> = lands under transportation and communication infrastructure; *ATRL*<sup>11</sup> = wetlands; *ATRL*<sup>12</sup> = disturbed lands; *ATRL*<sup>13</sup> = barren; bold denotes a strong correlation, *CATRli* > *Ccv* (0.5904 for type II territories). Source: Authors' development.

The yellow belt included three types of territories, namely, northern territories, Siberia, and the Far East, occupying over half of the territory of Russia, but only representing 12.3% of its agricultural land, where the land use was primarily rangeland. The variations in the acreage of rangelands strongly correlated with those of woodlands, other forest ranges, and wetlands (Table 5). The northern locus included the territories of Russia's northwest, the Ural region, and central Russia (i.e., north of Moscow). In these highly industrialized but less populated territories, we revealed strong correlations between the proportions of croplands and barren land, as well as between those of perennial plantings and disturbed lands. In the south, the yellow belt included Krasnodar, the principal breadbasket territory of Russia. The share of cropland in the composition of Krasnodar's land fund was 52.8%. Krasnodar is also one of Russia's most densely populated regions and is the most popular resort area. The analysis demonstrated high correlations between the proportions of cropland and perennial plantings, on one side, and the shares of residential and industrial lands and lands under transportation and communication infrastructure on the other.


**Table 5.** Correlation matrix for type III territories.

Note: *ATRLi* = centered log-ratio-transformed data: *ATRL*<sup>1</sup> = cropland; *ATRL*<sup>2</sup> = fallow; *ATRL*<sup>3</sup> = perennial plantings; *ATRL*<sup>4</sup> = hayfields; *ATRL*<sup>5</sup> = rangeland; *ATRL*<sup>6</sup> = woodlands; *ATRL*<sup>7</sup> = forest range; *ATRL*<sup>8</sup> = water reserve lands; *ATRL*<sup>9</sup> = residential and industrial lands; *ATRL*<sup>10</sup> = lands under transportation and communication infrastructure; *ATRL*<sup>11</sup> = wetlands; *ATRL*<sup>12</sup> = disturbed lands; *ATRL*<sup>13</sup> = barren; bold denotes a strong correlation, *CATRli* > *Ccv* (0.7458 for type III territories). Source: Authors' development.

Type IV comprised the territories with the lowest activity of agricultural lands. The scarcity of agricultural lands represented intercategory variations in the composition of the agricultural land fund. The strongest correlations were identified between various categories of agricultural lands, specifically, cropland and hayfields, on one side, and perennial plantings and rangeland on the other (Table 6). The composition of the agricultural land fund was also affected by the proportions of barren land (in Chukotka and Nenets), woodlands (in Leningrad and Murmansk), wetlands (in Murmansk), and water reserve lands (in Yamal-Nenets).



Note: *ATRLi* = centered log-ratio-transformed data: *ATRL*<sup>1</sup> = cropland; *ATRL*<sup>2</sup> = fallow; *ATRL*<sup>3</sup> = perennial plantings; *ATRL*<sup>4</sup> = hayfields; *ATRL*<sup>5</sup> = rangeland; *ATRL*<sup>6</sup> = woodlands; *ATRL*<sup>7</sup> = forest range; *ATRL*<sup>8</sup> = water reserve lands; *ATRL*<sup>9</sup> = residential and industrial lands; *ATRL*<sup>10</sup> = lands under transportation and communication infrastructure; *ATRL*<sup>11</sup> = wetlands; *ATRL*<sup>12</sup> = disturbed lands; *ATRL*<sup>13</sup> = barren; bold denotes a strong correlation, *CATRli* > *Ccv* (0.6293 for type IV territories). Source: Authors' development.

#### **4. Discussion**

The results, as expected, demonstrated that the compositions of the land funds in Russia vary across territories. Echoing Bichler et al. [100], Chu [101], Smith et al. [102], and Bakker et al. [103], we found that the distribution of agricultural lands is largely affected by natural factors, while agricultural lands are spread unevenly across the country. At a regional scale, belt-type concentrations of cropland suggest an agriculture-focused land distribution pattern in the southern and central areas of Russia. This is consistent with the observations of Rounsevell et al. [104] and White and Engelen [105,106], who revealed that agricultural land use tends to become concentrated in locations, reflecting the influence of natural factors and neighboring land distribution patterns. Nevertheless, in particular territories, the proportion of agricultural lands in the land funds do not match the type of agricultural land activity.

Emulating earlier studies by Mazurkin and Mihailova [70], Shishkina et al. [74], Mazurkin [69], and Buckett [71], we revealed that the application of a land activity parameter could result in creating a picture of land distribution patterns that are different from that which might be expected from the knowledge of the proportions of individual land categories. Therefore, land distribution change maps are not sufficient to capture specific finer-scale variations in the compositions of land funds at a regional scale. In Russia, land tenure and demand for land have been the principal economic proxies to map agricultural land distribution. According to Shagaida [107], the demand for agricultural land varies significantly across Russia's territories, depending on the degree of land consolidation. In the course of land reform, the previously dominant state farms have transformed the organizational form of their land use but still have persisted as the backbone of the agricultural sector [34,108]. In the embryonic land market in the 1990–2000s, the establishment of new land tenure patterns had not involved immediate changes in the distribution of land from big ex-Soviet agricultural enterprises to individual owners [107]. Since land certificates do not specify land plots, most of the shareowners have not withdrawn their land property from joint use by former collective farms. Over 70% of land in Russia is still used by large enterprises for rent, 25% is contributed to the capital of large enterprises, and only 4% is retained by private owners [109]. In the breadbasket southern and central European territories of Russia, large agricultural holding companies have aggregated even more agricultural land property when compared to the Soviet period [110].

To a large extent, the existing demand-based distribution matches the land activity map (Figure 3), as the highest demand for land is identified in the central parts of the country close to Moscow. This demand primarily exists due to non-agricultural businesses. For type I and II territories, this correlates well with the finding of strong links between the proportions of agricultural land categories, on one side, and those of residential, industrial, and infrastructure lands on the other. In type III southern locus (Krasnodar), Lerman and Shagaida [20] reported high demand for land among corporate farms. In that classification, type I and II territories are considered as less developed areas in terms of agricultural production (sometimes even as "agriculturally depressed regions" ([20] p. 20)), where corporate farms tend to reduce their holdings and abandon land plots. Our results, on the contrary, demonstrated that in the south of European Russia, where the concentration of croplands is the highest, agricultural land activity is lower compared to many other territories of the country.

In the territories where a high proportion of croplands coexist with low agricultural land activity, many of the variations in the composition of a land fund could be explained by socio-economic factors. Van de Steeg et al. [111] and Gärtner et al. [112] confirmed that the distribution of agricultural land strongly correlates with the level of rural development, proximity to economic and market centers, urbanization, and the demand for agricultural land from non-agricultural industries. Our study revealed correlations between the proportions of agricultural and urban lands across type I–III territories, which could represent losing agricultural land due to urban development. In type II territories, the compositions of agricultural land funds are more affected by urban development than the compositions of type I and III. These results supported the findings of Daniels [113], Su et al. [114], Yeh and Huang [115], and Dredge [116], i.e., the proximity to urban development can be a powerful predictor of changes in agricultural land use. Many scholars, including Parsipour et al. [117], Li et al. [118],

and Al-Kofahi et al. [119], among others, agree that the accelerating urbanization has been causing increasingly harmful effects on agricultural lands. In the case of Russia, we did not reveal the acceleration of agriculture land loss in urbanized type I–III territories. What was revealed, however, was the strengthening of the correlation between the variations in the compositions of agricultural land funds and residential, industrial, and infrastructure lands. As Zubair et al. [120] and Lucero and Tarlock [121] forecasted, such stronger associations would continue to put increasing pressure on agricultural lands and result in more fragmented agricultural land use in the future.

Along with urbanization, an orientation of a land fund composition towards agricultural production is determined by the population density [111,122]. In urbanized type I and II territories, agricultural land use is affected by the variations in the acreage of residential lands. In agricultureoriented Krasnodar, Rostov, and Stavropol, the changes in agricultural land fund compositions are mainly linked with those of lands for transportation and communication. This result was consistent with what Ramadani and Bytyqi [123], Li et al. [118], and Al-Kofahi et al. [119] reported when assessing the effects of more significant concentrations of the population on the lower proportions of agricultural lands in a land fund.

Reversely, Meyfroidt et al. [124] and Nguyen et al. [125] revealed that in the industrialized areas in Russia, where the density of population is lower, the concentration of abandoned lands is higher. There is an array of studies that have reported a link between industrial growth, changes in agricultural land distribution, and the degradation of farming opportunities internationally. Explicitly, Oyebanji et al. [126] confirmed the existence of a positive long-term relationship between industrialization and land loss in Nigeria. Deng and Li [127] revealed that the soil sealing effect has resulted from industrial and infrastructure construction in China, while Müller and Sikor [128], Milanova et al. [33], and Müller et al. [129] linked changes in agricultural land distribution and agricultural abandonment in EU countries with unfavorable environmental conditions due to increasing industrialization. The expansion of urban and industrial infrastructure not only triggers agriculture-to-urban and agriculture-to-industry land transfers but also leads to the overexploitation and degradation of remaining agricultural lands [127]. Many areas in Russia may soon face a reduction in farming opportunities due to various kinds of environmental pollution. Many experts tend to explain the unprecedented increase of barren land in Russia (by four million ha during the past two decades) by the intensive exploitation of mineral resources and industrial construction [39,130]. Kashtanov [131] and Dobrovolski [132] associated the expansion of industrial infrastructure with long-term and irreversible losses of cropland in Russia. In support of the earlier findings of Sorokin et al. [130] and Solgerel et al. [133] concerning the close relationships between industrial development and arable acreage, strong correlations between the proportions of croplands, perennial plantings, and industrial lands are revealed in both urbanized type I and II territories and sparsely populated yellow belt areas.

Distinct from urbanization, industrialization may affect agricultural land use in remote areas. According to Sorokin et al. [130], most of the abandoned lands are located in the north of Russia. This agrees well with our finding of strong correlations between the variations in the acreage of croplands, disturbed lands, and barren lands in the north locus of the yellow belt. Prishchepov et al. [39] and MacDonald et al. [134] also reported abandoned agricultural land concentrated in remote and isolated industrialized areas in northern Russia. Nakvasina et al. [135] claimed that the proximity to urban areas might be used as a critical criterion to transfer disturbed and barren lands back into agricultural use. However, we did not identify strong correlations between the variations in agricultural land fund compositions and residential lands for type III territories.

In diverse land activity patterns across the Russian territories, changes in the compositions of agricultural and non-agricultural land funds depend on the degree of industrial development. As mentioned by Postek et al. [136] and Prishchepov et al. [39], agricultural land loss due to increasing industrialization causes the fragmentation of arable lands as smaller locations with lower productivity. However, according to Popov [137], fragmentation is not a problem in agriculture-oriented areas due to the excessive lease of agricultural land. The issue is particularly topical in territories where

arable land is scarce, however [138,139]. Nefedova [140] reported that in northern and eastern parts of Russia, agricultural land distribution is extremely fragmented. Our results demonstrated that in the Russian North and Far East, low activity of cropland is coupled with the prevalence of hayfields in the composition of the agricultural land funds there. High intragroup correlations between the proportions of cropland, rangeland, hayfields, and perennial plantings in type IV territories confirm the observations of King and Burton [141], Tan et al. [142], and Dhakal and Khanal [143], i.e., the fragmentation results in the competition between the categories of agricultural lands.

We performed our analysis in the short-term, but it is commonly known that land transformations (particularly, for croplands and annual crops) can be rapid, whereas transformations are slower in grassland-livestock oriented areas and permanent crop areas. Nationally, the ongoing loss of croplands may not have an immediate effect on the agricultural output of Russia. Still, this represents enormous environmental, economic, and social costs that will be hard to absorb in terms of a long-term perspective [144]. Griewald et al. [34] and Hunt et al. [145] outlined five principal drivers of long-term change in agricultural land use in Russia, environmental drivers being one of them. Our findings would allow one to expect that the evolution of land-use change will be affected by the pressure exerted on ecosystems by various land management types [34]. While some authors, including Diputra and Baek [146] and Mahcene et al. [147], reported little evidence that industrialization causes a significant increase in disturbed land acreage, our results suggested that lower activity of agricultural land categories is correlated with a higher activity of barren land, disturbed land, and industrial land. Weaker, but still significant, correlations between the proportions of agricultural and industrial land categories are revealed in type I and II territories here. In type IV territories, the contributions of croplands and perennial plantings to regional land funds are also linked with variations in the acreage of barren lands.

Among the drivers of land-use change, in the long run, there are also economic, social, technological, and policy-related factors. Bukvareva et al. [148] stated that current land-use policies in Russia pay little attention to the environmental costs associated with the re-use of abandoned lands. In light of the economic recession that Russia has been experiencing since the mid-2010s, farmers tend to reinforce the exploitation of all available lands to ensure sufficient income inflow. Often, this is done regardless of whether some lands are of high environmental value or are socioeconomically marginal [29]. In the short-term perspective, we did not reveal an increase in the acreage of croplands due to the use of other categories of agricultural land. To some extent, however, the correlations between the proportions of agricultural land categories are identified in type III territories. In these yellow belt areas, land reclamation programs will require substantial investments for clearing forested land, liming, and other works. In the short-term, high reclamation costs along with poor soil quality may reduce expected economic returns [149]; however, in the long run, the incentives for reclamation may grow as both the availability and quality of croplands in type I and II territories degrade. Such a perspective highlights the need for a deeper investigation of the variations in land fund compositions within a sustainable agricultural land management approach as a component of the broader economic and environmental system [150,151].

#### **5. Conclusions**

In recent decades, there has been increasing concern for ensuring the effective utilization of agricultural land due to the limited area of highly productive arable land and the growing demand for food and farming products internationally. In Russia, an orientation of state policy towards the growth of agricultural production, along with a low level of environmental awareness among farmers, has impeded the prospects of sustainable land management as an integral aspect of land use planning. The degradation of agricultural lands due to irrational use has posed environmental, economic, and social threats to the national development objectives of land management in many territories of the country. As most studies in Russia have focused on land changes between the categories of agricultural land, the influence of agriculture-to-urban and agriculture-to-industry transfers has been downplayed.

We conducted this work, intending to study such variations by revealing the interdependencies between the proportions of agricultural land categories, on the one hand, and urban, industrial, and other types of land on the other. First, land distribution was mapped based on a share of agricultural lands in a composition of a land fund and, second, by a "land activity rating" of Russia's territories. Such a two-step approach to mapping allowed us to find that the proportions of agricultural lands in the composition of a land fund do not appropriately reveal the variations in the activities of agricultural land categories. In the territories, where agricultural lands dominated in the structure of a land fund, the agricultural land activity could be depressed by high proportions of non-agricultural lands. In urbanized and densely populated territories, the composition of the agricultural land fund was predominantly affected by the changes in the acreage of residential and industrial lands, as well as the lands for transportation and communication. In industrialized but underpopulated territories, the acreages of croplands and perennial plantings were strongly correlated with those of disturbed and barren lands. We also found that lower land activity tended to increase the variations within the agricultural land fund, which might indicate intercategory competition for more fertile, more productive, and better-located agricultural lands.

By establishing and testing the five-stage algorithm, we attempted to solve the scientific problem of low awareness in the causality between land-use processes and the composition of the land funds at regional scales. As distinguished from previous studies in the area, we investigated variations in the compositions of a land fund as interactions between the proportions of agricultural and non-agricultural lands. Practically, in territory-scale studies, such an approach might complement regionally adapted monitoring networks by targeting the mismatches between the cadaster-based mappings of agricultural land distributions and ranking-based activities of agricultural lands. Theoretically, such an algorithm allows one to capture the complex relationships of a variety of land categories and the resulting correlations between their proportions, therefore, being applicable for studying territorial land-use patterns, the simulation of agricultural land distribution systems, and the extrapolation of current trends into the future. Potentially, the algorithm is suitable for numerous locations. However, one of the limitations of the current study was that it used the Russian system of land statistics, which is built on thirteen land categories. Due to the different sources of land use data in different countries, an adjustment of the array of land categories to a national land reporting system is needed when implementing the method in a broader international context. Another limitation that could potentially challenge cross-country comparisons is the different sizes of territorial units. Russia's case demonstrates that this problem may arise even within one country, where territories substantially differ in size. In an attempt to overcome a data discrepancy obstacle, we conversed cadastral classification data into land-rating values. To address the diversity of territories, we used an agricultural land activity parameter. This allowed us to adjust agricultural and non-agricultural-oriented ranking systems to make them comparable. Nevertheless, further research is needed to assess to what extent the approach would be able to appropriately picture variations in agricultural land activity patterns in the conditions of information asymmetries among countries.

**Author Contributions:** V.E. conceptualized the research framework and wrote the paper; T.G. analyzed the data; A.I. performed the data collection. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research and the APC were funded by the International and Regional Studies Fund of the Ministry of Education of the People's Republic of China (grant number 19GBQY117).

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

#### **Appendix A**

**Table A1.** Land acreage data of the Central Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>1</sup> = Belgorod; *T*<sup>2</sup> = Bryansk; *T*<sup>3</sup> = Vladimir; *T*<sup>4</sup> = Voronezh; *T*<sup>5</sup> = Ivanovo; *T*<sup>6</sup> = Kaluga; *T*<sup>7</sup> = Kostroma; *T*<sup>8</sup> = Kursk; *T*<sup>9</sup> = Lipetsk; *T*<sup>10</sup> = Moscow Oblast; *T*<sup>11</sup> = Orel; *T*<sup>12</sup> = Ryazan; *T*<sup>13</sup> = Smolensk; *T*<sup>14</sup> = Tambov; *T*<sup>15</sup> = Tver; *T*<sup>16</sup> = Tula; *T*<sup>17</sup> = Yaroslavl; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A2.** Land acreage data of Northwestern Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>18</sup> = Karelia; *T*<sup>19</sup> = Komi; *T*<sup>20</sup> = Arkhangelsk; *T*<sup>21</sup> = Vologda; *T*<sup>22</sup> = Kaliningrad; *T*<sup>23</sup> = Leningrad; *T*<sup>24</sup> = Murmansk; *T*<sup>25</sup> = Novgorod; *T*<sup>26</sup> = Pskov; *T*<sup>27</sup> = Nenets; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A3.** Land acreage data of the Southern Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>28</sup> = Adygeya; *T*<sup>29</sup> = Kalmykia; *T*<sup>30</sup> = Crimea; *T*<sup>31</sup> = Krasnodar; *T*<sup>32</sup> = Astrakhan; *T*<sup>33</sup> = Volgograd; *T*<sup>34</sup> = Rostov; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A4.** Land acreage data of the North Caucasian Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>35</sup> = Dagestan; *T*<sup>36</sup> = Ingushetia; *T*<sup>37</sup> = Kabardino-Balkaria; *T*<sup>38</sup> = Karachaevo-Cherkessia; *T*<sup>39</sup> = North Osetia-Alania; *T*<sup>40</sup> = Chechnya; *T*<sup>41</sup> = Stavropol; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A5.** Land acreage data of the Volga Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>42</sup> = Bashkortostan; *T*<sup>43</sup> = Mari El; *T*<sup>44</sup> = Mordovia; *T*<sup>45</sup> = Tatarstan; *T*<sup>46</sup> = Udmurtia; *T*<sup>47</sup> = Chuvashia; *T*<sup>48</sup> = Perm; *T*<sup>49</sup> = Kirov; *T*<sup>50</sup> = Nizhny Novgorod; *T*<sup>51</sup> = Orenburg; *T*<sup>52</sup> = Penza; *T*<sup>53</sup> = Samara; *T*<sup>54</sup> = Saratov; *T*<sup>55</sup> = Ulyanovsk; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A6.** Land acreage data of the Ural Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>56</sup> = Kurgan; *T*<sup>57</sup> = Sverdlovsk; *T*<sup>58</sup> = Tyumen; *T*<sup>59</sup> = Chelyabinsk; *T*<sup>60</sup> = Khanty-Mansi; *T*<sup>61</sup> = Yamal-Nenets; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.


**Table A7.** Land acreage data of the Siberian Federal District in thousand hectares. Mean values for 2010–2018.

Note: *T*<sup>62</sup> = Altay Republic; *T*<sup>63</sup> = Buryatia; *T*<sup>64</sup> = Tyva; *T*<sup>65</sup> = Khakasia; *T*<sup>66</sup> = Altay; *T*<sup>67</sup> = Zabaikalsk; *T*<sup>68</sup> = Krasnoyarsk; *T*<sup>69</sup> = Irkutsk; *T*<sup>70</sup> = Kemerovo; *T*<sup>71</sup> = Novosibirsk; *T*<sup>72</sup> = Omsk; *T*<sup>73</sup> = Tomsk; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.

**Table A8.** Land acreage data of the Far Eastern Federal District in thousand hectares. Mean values for 2010–2018.


Note: *T*<sup>74</sup> = Sakha Yakutia; *T*<sup>75</sup> = Kamchatka; *T*<sup>76</sup> = Primorye; *T*<sup>77</sup> = Khabarovsk; *T*<sup>78</sup> = Amur; *T*<sup>79</sup> = Magadan; *T*<sup>80</sup> = Sakhalin; *T*<sup>81</sup> = Jewish AO; *T*<sup>82</sup> = Chukotka; *L*(1–13) = acreage of *Li* category, thousand hectares: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren. Source: Authors' development.


**Table A9.** Activity per land category in Russia, Central Federal District. Mean values for 2010–2018.

transportation

change or insignificant change. Source: Authors' development.

 and

communication

infrastructure; *L*11 =

wetlands; *L*12 =

disturbed lands; *L*13 = barren; *VL*(1–13) =

variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no


**Table A10.** Activity per land category in Russia, Northwestern Federal District. Mean values for 2010–2018.

Note: *T*<sup>18</sup> = Karelia; *T*<sup>19</sup> = Komi; *T*<sup>20</sup> = Arkhangelsk; *T*<sup>21</sup> = Vologda; *T*<sup>22</sup> = Kaliningrad; *T*<sup>23</sup> = Leningrad; *T*<sup>24</sup> = Murmansk; *T*<sup>25</sup> = Novgorod; *T*<sup>26</sup> = Pskov; *T*<sup>27</sup> = Nenets; *L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.

**Table A11.** Activity per land category in Russia, Southern Federal District. Mean values for 2010–2018.



**Table A11.** *Cont.*

Note: *T*<sup>28</sup> = Adygeya; *T*<sup>29</sup> = Kalmykia; *T*<sup>30</sup> = Crimea; *T*<sup>31</sup> = Krasnodar; *T*<sup>32</sup> = Astrakhan; *T*<sup>33</sup> = Volgograd; *T*<sup>34</sup> = Rostov; *L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.

**Table A12.** Activity per land category in Russia, North Caucasian Federal District. Mean values for 2010–2018.


Note: *T*<sup>35</sup> = Dagestan; *T*<sup>36</sup> = Ingushetia; *T*<sup>37</sup> = Kabardino-Balkaria; *T*<sup>38</sup> = Karachaevo-Cherkessia; *T*<sup>39</sup> = North Osetia-Alania; *T*<sup>40</sup> = Chechnya; *T*<sup>41</sup> = Stavropol; *L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.



*L*3 =

perennial plantings; *L*4 =

transportation

change or insignificant change. Source: Authors' development.

 and

communication

infrastructure; *L*11 =

wetlands; *L*12 =

disturbed lands; *L*13 = barren; *VL*(1–13) =

hayfields; *L*5 =

rangeland; *L*6 =

woodlands; *L*7 = forest range; *L*8 = water reserve lands; *L*9 =

residential and industrial lands; *L*10 = lands under

variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no


**Table A14.** Activity per land category in Russia, Ural Federal District. Mean values for 2010–2018.

Note: *T*<sup>56</sup> = Kurgan; *T*<sup>57</sup> = Sverdlovsk; *T*<sup>58</sup> = Tyumen; *T*<sup>59</sup> = Chelyabinsk; *T*<sup>60</sup> = Khanty-Mansi; *T*<sup>61</sup> = Yamal-Nenets; *L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.

**Table A15.** Activity per land category in Russia, Siberian Federal District. Mean values for 2010–2018.


**Table A15.** *Cont.*


Note: *T*<sup>62</sup> = Altay Republic; *T*<sup>63</sup> = Buryatia; *T*<sup>64</sup> = Tyva; *T*<sup>65</sup> = Khakasia; *T*<sup>66</sup> = Altay; *T*<sup>67</sup> = Zabaikalsk; *T*<sup>68</sup> = Krasnoyarsk; *T*<sup>69</sup> = Irkutsk; *T*<sup>70</sup> = Kemerovo; *T*<sup>71</sup> = Novosibirsk; *T*<sup>72</sup> = Omsk; *T*<sup>73</sup> = Tomsk;*L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.

**Table A16.** Activity per land category in Russia, Far Eastern Federal District. Mean values for 2010–2018.


Note: *T*<sup>74</sup> = Sakha Yakutia; *T*<sup>75</sup> = Kamchatka; *T*<sup>76</sup> = Primorye; *T*<sup>77</sup> = Khabarovsk; *T*<sup>78</sup> = Amur; *T*<sup>79</sup> = Magadan; *T*<sup>80</sup> = Sakhalin; *T*<sup>81</sup> = Jewish AO; *T*<sup>82</sup> = Chukotka; *L*(1–13) = portion of *Li* category in a composition of the land fund in *Tj* territory, percentage: *L*<sup>1</sup> = cropland; *L*<sup>2</sup> = fallow; *L*<sup>3</sup> = perennial plantings; *L*<sup>4</sup> = hayfields; *L*<sup>5</sup> = rangeland; *L*<sup>6</sup> = woodlands; *L*<sup>7</sup> = forest range; *L*<sup>8</sup> = water reserve lands; *L*<sup>9</sup> = residential and industrial lands; *L*<sup>10</sup> = lands under transportation and communication infrastructure; *L*<sup>11</sup> = wetlands; *L*<sup>12</sup> = disturbed lands; *L*<sup>13</sup> = barren; *VL*(1–13) = variability of *L*(1–13), i.e., change in 2018 compared to 2010; "-" = no change or insignificant change. Source: Authors' development.

#### **Appendix C**


**Table A17.** Ranking of *Tj* territories on land activity, Central Federal District.

Note: *T*<sup>1</sup> = Belgorod; *T*<sup>2</sup> = Bryansk; *T*<sup>3</sup> = Vladimir; *T*<sup>4</sup> = Voronezh; *T*<sup>5</sup> = Ivanovo; *T*<sup>6</sup> = Kaluga; *T*<sup>7</sup> = Kostroma; *T*<sup>8</sup> = Kursk; *T*<sup>9</sup> = Lipetsk; *T*<sup>10</sup> = Moscow Oblast; *T*<sup>11</sup> = Orel; *T*<sup>12</sup> = Ryazan; *T*<sup>13</sup> = Smolensk; *T*<sup>14</sup> = Tambov; *T*<sup>15</sup> = Tver; *T*<sup>16</sup> = Tula; *T*<sup>17</sup> = Yaroslavl; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.

**Table A18.** Ranking of *Tj* territories on land activity, Northwestern Federal District.


Note: *T*<sup>18</sup> = Karelia; *T*<sup>19</sup> = Komi; *T*<sup>20</sup> = Arkhangelsk; *T*<sup>21</sup> = Vologda; *T*<sup>22</sup> = Kaliningrad; *T*<sup>23</sup> = Leningrad; *T*<sup>24</sup> = Murmansk; *T*<sup>25</sup> = Novgorod; *T*<sup>26</sup> = Pskov; *T*<sup>27</sup> = Nenets; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.


**Table A19.** Ranking of *Tj* territories on land activity, Southern Federal District.

Note: *T*<sup>28</sup> = Adygeya; *T*<sup>29</sup> = Kalmykia; *T*<sup>30</sup> = Crimea; *T*<sup>31</sup> = Krasnodar; *T*<sup>32</sup> = Astrakhan; *T*<sup>33</sup> = Volgograd; *T*<sup>34</sup> = Rostov; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.

**Table A20.** Ranking of *Tj* territories on land activity, North Caucasian Federal District.


Note: *T*<sup>35</sup> = Dagestan; *T*<sup>36</sup> = Ingushetia; *T*<sup>37</sup> = Kabardino-Balkaria; *T*<sup>38</sup> = Karachaevo-Cherkessia; *T*<sup>39</sup> = North Osetia-Alania; *T*<sup>40</sup> = Chechnya; *T*<sup>41</sup> = Stavropol; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.

**Table A21.** Ranking of *Tj* territories on land activity, Volga Federal District.


Note: *T*<sup>42</sup> = Bashkortostan; *T*<sup>43</sup> = Mari El; *T*<sup>44</sup> = Mordovia; *T*<sup>45</sup> = Tatarstan; *T*<sup>46</sup> = Udmurtia; *T*<sup>47</sup> = Chuvashia; *T*<sup>48</sup> = Perm; *T*<sup>49</sup> = Kirov; *T*<sup>50</sup> = Nizhny Novgorod; *T*<sup>51</sup> = Orenburg; *T*<sup>52</sup> = Penza; *T*<sup>53</sup> = Samara; *T*<sup>54</sup> = Saratov; *T*<sup>55</sup> = Ulyanovsk; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.


**Table A22.** Ranking of *Tj* territories on land activity, Ural Federal District.

Note: *T*<sup>56</sup> = Kurgan; *T*<sup>57</sup> = Sverdlovsk; *T*<sup>58</sup> = Tyumen; *T*<sup>59</sup> = Chelyabinsk; *T*<sup>60</sup> = Khanty-Mansi; *T*<sup>61</sup> = Yamal-Nenets; *R*(1–13) =ranks of land activity per land categories: *R*<sup>1</sup> =cropland; *R*<sup>2</sup> =fallow; *R*<sup>3</sup> =perennial plantings; *R*<sup>4</sup> =hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.

**Table A23.** Ranking of *Tj* territories on land activity, Siberian Federal District.


Note: *T*<sup>62</sup> = Altay Republic; *T*<sup>63</sup> = Buryatia; *T*<sup>64</sup> = Tyva; *T*<sup>65</sup> = Khakasia; *T*<sup>66</sup> = Altay; *T*<sup>67</sup> = Zabaikalsk; *T*<sup>68</sup> = Krasnoyarsk; *T*<sup>69</sup> = Irkutsk; *T*<sup>70</sup> = Kemerovo; *T*<sup>71</sup> = Novosibirsk; *T*<sup>72</sup> = Omsk; *T*<sup>73</sup> = Tomsk; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.

**Table A24.** Ranking of *Tj* territories on land activity, Far Eastern Federal District.


Note: *T*<sup>74</sup> = Sakha Yakutia; *T*<sup>75</sup> = Kamchatka; *T*<sup>76</sup> = Primorye; *T*<sup>77</sup> = Khabarovsk; *T*<sup>78</sup> = Amur; *T*<sup>79</sup> = Magadan; *T*<sup>80</sup> = Sakhalin; *T*<sup>81</sup> = Jewish AO; *T*82= Chukotka; *R*(1–13) = ranks of land activity per land categories: *R*<sup>1</sup> = cropland; *R*<sup>2</sup> = fallow; *R*<sup>3</sup> = perennial plantings; *R*<sup>4</sup> = hayfields; *R*<sup>5</sup> = rangeland; *R*<sup>6</sup> = woodlands; *R*<sup>7</sup> = forest range; *R*<sup>8</sup> = water reserve lands; *R*<sup>9</sup> = residential and industrial lands; *R*<sup>10</sup> = lands under transportation and communication infrastructure; *R*<sup>11</sup> = wetlands; *R*<sup>12</sup> = disturbed lands; *R*<sup>13</sup> = barren. Source: Authors' development.



#### *Land* **2020**, *9*, 201

Rj9



territories per districts; *Rji* is averaged in respect to individual rankings *Rji* in *Tj* territories per districts;


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Land* Editorial Office E-mail: land@mdpi.com www.mdpi.com/journal/land

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-4050-4