*Article* **Water Regulation Ecosystem Services of Multifunctional Landscape Dominated by Monoculture Plantations**

**Yudha Kristanto 1,2,\*, Suria Tarigan 3, Tania June 4, Enni Dwi Wahjunie <sup>3</sup> and Bambang Sulistyantara <sup>5</sup>**


**Abstract:** Meeting the growing demand for agricultural production while preserving water regulation ecosystem services (WRES) is a challenge. One way to preserve WRES is by adopting multifunctional landscape approach. Hence, the main objective was to evaluate the role of forest patches (FP) in preserving WRES in tropical landscapes dominated by oil palm plantations. The SWAT model was used to evaluate the essential WRES, such as water yield (WYLD), soil water (SW), surface runoff (SURQ), groundwater recharge (GWR), and evapotranspiration (AET). Due to a compaction, soils in monoculture plantation have higher bulk density and lower porosity and water retention, which decrease WRES. Conserving FP among oil palms evidently improves WRES, such as decreasing SURQ and rain season WYLD and increasing GWR, SW, AET, and dry season WLYD. FP has sponge-like properties by storing water to increase water availability, and pump-like properties by evaporating water to stabilize the microclimate. Mature oil palm also has pump-like properties to maintain productivity. However, it does not have sponge-like properties that make water use more significant than the stored water. Consequently, a multifunctional landscape could enhance WRES of forest patches and synergize it with provisioning ecosystem services of oil palm plantations.

**Keywords:** evapotranspiration; groundwater recharge; soil compaction; surface runoff; soil water retention

#### **1. Introduction**

Local, countrywide, and global consumption of processed palm oil products, such as food, bioenergy, and oleochemicals, is growing along with population growth. This phenomenon implies changing the tropical landscapes into oil palm plantations (*Elaeis guineensis* Jacq.) that are intensively cultivated in monocultures [1]. The area of oil palm plantations has rapidly grown and is expected to increase in the coming years [2,3]. On the one hand, the development of oil palm plantations rapidly contributes to the local, regional, and national economies [4]. In addition, oil palm plantations also positively impact increasing access to the basic needs of local communities, such as school facilities, health facilities, road networks, and electricity networks. However, if the development of oil palm plantations is not appropriately managed, it can cause the degradation of ecosystem services, and reduce the ecological function of a landscape and harms the socio-economic community [5].

Previous research has examined the transformation of natural ecosystem into agroecosystem, especially monoculture plantations on ecosystem services changes [2]. One of the

**Citation:** Kristanto, Y.; Tarigan, S.; June, T.; Wahjunie, E.D.; Sulistyantara, B. Water Regulation Ecosystem Services of Multifunctional Landscape Dominated by Monoculture Plantations. *Land* **2022**, *11*, 818. https://doi.org/10.3390/ land11060818

Academic Editors: Chiara Piccini, Rosa Francaviglia and Richard Cruse

Received: 2 April 2022 Accepted: 28 May 2022 Published: 31 May 2022

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

**Copyright:** © 2022 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/).

noticeable changes in ecosystem services in agroecosystem is the change in water regulation ecosystem services (WRES) caused by soil compaction due to intensive land management. These phenomena encompass a decrease in infiltration [6], an increase in surface runoff in the rainy season [3], and a decrease in groundwater availability in the dry season [5]. Therefore, meeting the growing demand for oil palm production while maintaining WRES is a challenge faced by a landscape of oil palm plantations related to sustainable and climate smart agriculture. Ecosystem services support human well-being through provision, regulation, and culture formed by natural and manufactured ecosystem structures and processes. WRES are essential ecosystem services in sustainable development, referring to the quantity, quality, and time of water stored in and out of an ecosystem [7]. WRES benefits are freshwater supply, flood and drought protection, electrical power generation, irrigation, and aquatic ecosystem maintenance.

One way to balance landscape WRES is by applying a multifunctional landscape approach through preserving the remnants of forest patches among oil palm plantations. These forest patches can be in secondary dryland forests, riparian vegetation, or agroforestry systems designated as high conservation areas (HCV). A *multifunctional landscape* is a landscape that can serve multiple ecosystem services simultaneously, not only for provisioning services, but also for regulation, cultural, and support services [8]. The multifunctional landscape is a more realistic soil and water conservation approach to optimize ecosystem services, where forest patches serve water regulation and the other regulation services, while oil palm plantations serve crop production. Furthermore, the multifunctional landscape maps the potential trade-offs and win-win synergies between each ecosystem service, especially provisioning services, which often conflict with regulation, cultural, and support services [9].

Since forest patches are essential to maintain the balance of landscape-scale WRES, it is necessary to model and evaluate WRES on oil-palm-dominated landscapes assisted by dynamic models, such as the Soil and Water Assessment Tools (SWAT). The strength of SWAT in WRES simulation is that SWAT has complex parameters and model structures, which can simulate the influence of physical processes in the soil–plant–atmosphere continuum that affects WRES on the watershed scale. However, field observation cannot measure several SWAT parameters, so the uncertainty that arises from parameter justification needs to be evaluated and considered in selecting the output model to be interpreted. In addition, previous research on oil palm ecosystem services generally does not differentiate between mature and young oil palms, especially in the same tropical landscape [3,6]. Oil palm growths affect the WRES changes through soil compaction due to intensive cultivation activities and canopy development that affect hydrological parameters at the hydrologic response unit (HRU) scale.

Because of soil compaction, WRES evaluation in this study needs to consider the soil hydrological characteristics and soil water retention curve (SWRC) at the smallest analysis unit, which significantly affect the WRES due to different land use and crop dynamics. The relationship between SWRC and WRES occurs at the HRU scale, where the information from SWRC is inputted into the .sol database to update the HRU value. SWRC was closely related to calculating the retention parameter of curve number and the range of soil moisture dynamics that affect surface runoff, evapotranspiration, and soil water storage. When the HRU as the smallest unit of analysis is updated, the subbasin parameters are automatically updated through the routing mechanism. Therefore, this study answers the main research question: how do forest patches preserve WRES in landscapes dominated by oil palm plantations? Consequently, two objectives can be drawn in this study: (i) analyze the soil water retention characteristics of each land use due to land management, and (ii) evaluate the role of forest patches among oil palm plantations in preserving landscape-scale WRES.

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

#### *2.1. Study Area*

The study area is a tropical lowland located in Jambi Province, Sumatra, Indonesia. This area has an equatorial rainfall type, characterized by two peaks of the rainy season (December and April) and one dry season (July) in one year. This area also has the same soil type, Hapludults, with a sandy clay loam texture (50% sand, 30% clay, and 20% silt) with small topography variations. The land use is dominated by oil palm plantations that are cultivated in monoculture, both mature plants (MOP) and young plants (YOP) (Figure 1). Among oil palm plantations, there are remnants of forest patches (FP) in the form of secondary dryland forest, riparian vegetation, and agroforestry plot separated from their primary ecosystem. The agroforestry plot (AGF ex-MOP) was planted in 2013 from former mature oil palm plantations.

**Figure 1.** Study area.

A micro-watershed with 13 sub-basins is delineated as a system boundary for WRES modeling. The micro-watershed has an area of 19.1 square kilometers, and lies between 1◦55 38.7–1◦58 48.3 S and 103◦11 50.1–103◦15 34.8 E, with an elevation range of 27–106 m above sea level. An automatic water level recorder (AWLR) is installed in the watershed outlet to support model simulation. The water level data are converted into streamflow through river morphometry measurements and derived into a rating curve equation using an exponential regression model. The rating curve equation obtained from the measurement is Q = 0.091 exp (2.309 H), where Q is the river discharge and H is the water level. In addition, an automatic weather station (AWS) is installed to record hourly and daily meteorological parameters. Meteorological measurements were carried out from January 2015 to January 2021 as model input, while river discharge measurements were carried out from November 2020 to January 2021 for model calibration and validation.

#### *2.2. Soil Water Retention Characteristics*

Undisturbed soil samples at a 0–15 cm depth were collected using soil ring samples with a 7.5 cm diameter and 5 cm height. The total sample for each land use was 8, consisting of 4 different locations and 2 replications for each location. There is a 20 × 20-m plot at each location divided into four quadrants with a size of 10 × 10 m. The soil was sampled in quadrants II and IV, which were diagonal to each other so that, at each location, there were two replications. Because the landscape tends to be homogeneous in terms of elevation, topography, and soil types, we assume that differences in ecosystems are only caused by

land use, so that differences in land use can represent different ecosystems at the landscape level. The undisturbed soil samples were used to determine soil hydrological properties, such as bulk density and porosity, and soil water retention characteristics at a particular suction matrix, such as pF 0 (saturated water content), pF 1, pF 2, pF 2.54 (field capacity), and pF 4.2 (permanent wilting point). Analysis of variance (ANOVA) at the 95% level was used to analyze differences in soil hydrological properties between land uses in the study area. The post hoc test with Duncan's multiple range test (DMRT) is conducted if the ANOVA showed a significant difference (*p*-Value < 0.05).

The observed soil hydrological properties are then used to model the soil water retention curve (SWRC). SWRC is a curve that describes the characteristics of soil water retention by defining the relationship between volumetric water content (θ) and matrix suction (ψ). SWRC is generally defined using a mathematical equation or pedotransfer function (PTF). The most used mathematical equation for SWRC modeling is the van Genuchten equation, as follows in Equation (1) [10]. The modeled SWRC was then calibrated using the observed moisture content at the pF 1, pF 2, pF 2.54, and pF 4.2.

$$\Theta(\mathbf{h}) = \Theta \mathbf{r} + \frac{\theta \mathbf{s} - \ \theta \mathbf{r}}{\left(1 + |\alpha \mathbf{H}|^{\mathbf{n}}\right)^{\mathbf{m}}} \tag{1}$$

where θ(h) is the soil moisture content (%*v*/*v*), θr is the residual water content (%*v*/*v*), θs is the saturated water content (value equal to the total porosity) (%*v*/*v*), α is the parameter related to the air entry value into the saturated soil (pF) (ψα = suction where the saturated soil goes through desaturation or air begins to enter the soil pores), H is the matrix suction in a logarithmic scale (pF), and n and m are parameters related to the slope of the curve at the inflection point (ψ > ψα). The slope of the curve (S) (1/pF) as a function of n and m can be calculated based on Equation (2) [11]:

$$\mathbf{S} = -\mathbf{n} \begin{pmatrix} \theta \mathbf{s} - \ \theta \mathbf{r} \end{pmatrix} \begin{bmatrix} 1 + \frac{1}{\mathbf{m}} \end{bmatrix}^{-(1+\mathbf{m})} \tag{2}$$

Furthermore, soil data obtained from sampling and laboratory analysis besides SWRC, such as soil permeability, texture, and soil organic matter, were used as inputs for the .sol database inside the SWAT model. The observed soil data are beneficial for WRES simulation, such as determining the initial abstraction, curve number calculation, soil moisture modeling, and evapotranspiration modeling.

#### *2.3. Evaluation of Water Regulation Ecosystem Services*

#### 2.3.1. Simulation of Water Regulation Ecosystem Services Using SWAT Model

The Soil and Water Assessment Tools (SWAT) is a physically based, computationally efficient, spatially semi-distributed, and temporally continuous watershed-scale ecosystem services model [12]. This model was developed to simulate various essential ecosystem services related to soil and water in the watershed system, such as water regulation, nutrient retention, and erosion prevention [13,14]. Compared to other WRES models, SWAT integrates the hydrological model with soil and land management attributes, such as irrigation, drainage, fertilization, tillage, and pesticides, and integrates the hydrological model with the crop growth model to simulate dynamic WRES related to crop growth phases [15]. SWAT also simulates various soil and water conservation options and ecological disturbance scenarios, such as land-use change and climate change (LUCCC), that affect ecosystem services [16]. It makes the SWAT results widely used as decision-making tools related to natural resources and environmental management, such as flood and drought mitigation, hydroelectric power generation, and LUCCC impact assessment [15].

SWAT simulates WRES, such as soil moisture (SW), surface runoff (SURQ), lateral flow (LAT), baseflow (BFO), and actual evapotranspiration (AET) from observed precipitation (PRECIP) data based on conservation of mass [14], as described by the Equation (3):

$$\text{SW}\_{\text{t}} = \text{SW}\_{0} + \sum\_{i=1}^{\text{t}} (\text{PRECIP} - \text{SURQ} - \text{AET} - \text{LAT} - \text{BFO}) \tag{3}$$

TheWRES simulation follows three stages: (i) preprocess, (ii) run model, and (iii) calibration, validation, and sensitivity analysis, and (iv) model uncertainty test.

(i) Preprocess

The preprocess stage includes watersheds, sub-basin, and river networks delineation from elevation data and defining HRU from land cover, soil type, and slope class. This stage resulted in a 19.1 km<sup>2</sup> study area enclosed within the micro-watershed boundary. SWAT divides a basin into sub-watersheds and divides a sub-watershed into hydrological response units (HRU) as the smallest unit of analysis. The defined HRUs are 105 HRUs, consisting of four land covers, one soil type, three slope classes, and thirteen subbasins combinations. HRU explains the spatial heterogeneity of WRES within the watershed and improves the accuracy of WRES modeling for any combination of land use, land management, vegetation type, soil, topography, and climate [17–19].

The spatial data needed to run the SWAT model is a digital elevation model (DEM) with a resolution of 8 m from the Indonesian Geospatial Information Agency, land use derived from Landsat-8 OLI with a resolution of 30 m, and soil map unit with a scale of 1:50,000 from Indonesian Center for Agricultural Land Resources Research and Development. The temporal data needed are daily meteorological data, including rainfall, air temperature, solar radiation, wind speed, and humidity. The preprocessing stage also includes attribute data input, which includes land cover (.mgt), soil physical properties (.sol), rainfall (.pcp), and potential evapotranspiration (.pet).


SWAT simulate WRES on each HRU and accumulate the WRES from HRU-scale to landscape-scale by a flow routing mechanism [20]. SWAT provides several methods for modeling WRES and flow routing, where users can choose which combination of methods best suits the characteristics of the study area and the availability of input data. One commonly used method for rainfall –runoff modeling related to WRES is the Soil Conservation Services-Curve Number (SCS-CN) [21]. SCS-CN is a powerful method of generating surface runoff, calculated based on Equation (4):

$$\text{SURQ} = \frac{\left(\text{PRECIP} - \text{I}\_{\text{a}}\right)^{2}}{\left(\text{PRECIP} - \text{I}\_{\text{a}} + \text{S}\right)^{\prime}}, \text{ where } \text{S} = 254 \left(\frac{100}{\text{CN}} - 1\right) \tag{4}$$

where SURQ is surface runoff (mm), PRECIP is precipitation (mm), S is retention parameter (mm), and Ia is an initial abstraction (mm). Initial abstraction is generally assumed 0.2 of the retention parameter, and SURQ only occurs when P > Ia (SURQ = 0 if P ≤ Ia). The three main processes considered in initial abstraction are rainfall interception, surface depression storage, and infiltration before the surface runoff. The retention parameter (S) can be approximated as the curve number (CN) function.

Due to land use, hydrologic soil group, and land management differences, the CN value varies spatially. The daily CN value will also vary temporally by considering the antecedent moisture condition (AMC): CN1—dry (wilting point), CN2—average moisture, and CN3—wet (saturated). The CN2 value for each land use, HSG, and land management was taken from the reference table [22], while CN1 and CN3 were calculated from CN2 [14]. SWAT provides two CN methods, original CN and modified or plants evapotranspiration curve number. In a previous study, [23] concluded that these two methods simulated streamflow with equally good performance but differed in simulating soil moisture and evapotranspiration in the lowland tropical landscape. The original CN was chosen in this study because it has a better performance in simulating various elements of WRES than CN-ET, including river discharge, soil water storage, and evapotranspiration.

#### • Evapotranspiration Modeling

Evapotranspiration as WRES also plays an important role in managing water resources, such as irrigation, soil–vegetation–atmosphere interactions, and spatial–temporal ecosystem productivity. The calculation of evapotranspiration by SWAT is based on the water continuity equation (Equation (3)) and potential evapotranspiration (PET). The PET model used to derive the actual evapotranspiration (AET) is the Penman–Monteith equation [14,24]. The PET calculation is automatically carried out by the SWAT model, which requires a database of maximum and minimum air temperature (.tmp), solar radiation (.slr), wind speed (.wnd), air humidity (.rhu), and crop parameters, such as leaf area index (LAI). The above meteorological datasets (.tmp, .slr, .wnd, and .rhu) were obtained from direct measurements in the field through weather monitoring with AWS. In addition, LAI data were also obtained from sampling using a hemisphere camera. After PET calculations, SWAT estimated AET as the sum of the canopy evaporation, soil evaporation, and plant transpiration [25]. Plant transpiration was calculated as a function of LAI, canopy evaporation was calculated as a function of rainfall interception, and soil evaporation was calculated as soil moisture [14].

#### (iii) Calibration, Validation, and Sensitivity Analysis

Simulations were carried out from 2015 to 2021, where 2015 was used to warm up the model. Calibration and validation are based on observed streamflow because they are easy to measure and cost-effective compared to soil moisture and evapotranspiration measurements [26]. The first half discharge data are used for the calibration, and the rest are used for validation. The calibrated parameters are parameters that cannot be measured directly in the field, such as groundwater and routing parameters. In contrast, parameters that can be measured directly, such as soil parameters and leaf area index are not calibrated. Compared with manual calibration, which takes a long time and fails to identify parameter sensitivity, this study uses automatic calibration based on the sequential uncertainty fitting-2 (SUFI-2) algorithm using SWAT-CUP software [27]. Sensitivity analysis was conducted to determine the response of changes in SWAT parameters to the significance of output changes and explore all possible combinations of model parameters to investigate output responses related to interactions between parameters [28]. The combination of model parameters and possible outputs are paired and sampled using the Latin hypercube sampling (LHS) to map their interactions and measure the output uncertainty caused by each parameter combination [27].

The Nash–Sutcliffe efficiency (NSE) is a selected model reliability indicator for the objective function during the parameter calibration. NSE value is used to measure how accurately the model's simulation results can describe the observation data. The NSE values range from -∞, which indicates that the model is highly inaccurate, to 1, which indicates that the model is highly accurate.

$$\text{NSE} = 1 - \left[ \frac{\sum\_{i=1}^{n} \left( \mathbf{Y}\_{\text{i}}^{\text{obs}} - \mathbf{Y}\_{\text{i}}^{\text{sim}} \right)^{2}}{\sum\_{i=1}^{n} \left( \mathbf{Y}\_{\text{i}}^{\text{obs}} - \overline{\mathbf{Y}^{\text{obs}}} \right)^{2}} \right] \tag{5}$$

where Yobs is observation data and Ysim is simulation data.

The model's reliability can also be evidenced by the coefficient of determination (R-squared) value. R-squared values range from 0, indicating that the model is highly inaccurate, to 1, which indicates that the model is very accurate. There are no absolute criteria for assessing the reliability of the hydrological model described in the literature. However, some criteria are commonly used, such as the NSE criteria by Moriasi [29] and R-squared criteria by Ayele [30].

#### (iv) Model Uncertainty Test

Evaluation of the model reliability is not enough to ensure that the SWAT outputs are genuinely accurate and interpreted directly. Furthermore, it is also necessary to evaluate the model uncertainty that arises due to the complexity of the SWAT structures and parameters justification. Therefore, the SUFI-2 algorithm in SWAT-CUP introduces statistical indicators to investigate the structural uncertainties associated with model simulations [27]. The uncertainty that arises during the parameter calibration is measured by the p-Factor, which is the percentage of observed data that is within 95% predictive uncertainty between the 2.5 and 97.5 percentiles (95PPU), and the r-Factor, which indicates the thickness of the mean of 95PPU divided by the standard deviation of the observed data [31]. Besides looking for high NSE and R-squared, it is also necessary to obtain the largest possible p-Factor and the smallest possible r-Factor. The uncertainty of the model is acceptable if the p-Factor > 0.7 and the r-Factor < 1.5 [27].

#### 2.3.2. Model Limitation

Parameter optimization during the calibration process can produce identical streamflow output with observational data regardless of how the best-fit parameters affect other WRES imprecision. However, because this research is related to the WRES assessment, the interpretation of the model is based not only on streamflow outputs, but also on other WRES, such as soil water storage and actual evapotranspiration. This study obtained precipitation as WRES input and other meteorological data for ETP calculation from the automatic weather station. Due to the limitations of time-series observations of soil moisture, we used soil hydrological properties and soil water retention curve (SWRC) observations from soil sampling and laboratory analysis. We linked the information from SWRC with the SWAT model by updating the .sol database for each HRU as soil moisture modeling inputs.

SWAT simulates soil moisture for each HRU as soil water storage (mm) in the range of available water content (AWC) between permanent wilting point (WP) and field capacity (FC). To obtain %*v*/*v* AWC, SWAT divides the soil water storage (mm) by soil depth (SOL\_Z) and adds this result with WP. Based on the information of FC, AWC, and WP from SWRC, the results of the soil moisture from the SWAT model are still within the AWC range following AWC observations on each land use. Finally, we consider the actual evapotranspiration as the "residual" component of the modeling based on the water balance equation (AET = PRECIP − Q − ΔSW). Therefore, the reliability of meteorological observation, SWRC observation, streamflow modeling, and SW modeling would affect the reliability of AET. If we could appropriately simulate the streamflow and soil moisture, then the AET value can also be relied upon in the future WRES evaluation.

#### 2.3.3. Water Regulation Ecosystem Services Indicators

One of the further challenges of evaluating WRES is determining the essential indicators based on the SWAT outputs. In general, precipitation is distributed into three flow elements: surface runoff (SURQ), groundwater recharge (GWR), and actual evapotranspiration (GWR), and the remainder is stored as soil water storage (SW).

$$\text{PRECIP}\_{\text{n}} = \text{SURQ}\_{\text{n}} + \text{GWR}\_{\text{n}} + \text{AET}\_{\text{n}} + \text{SW}\_{\text{n}} \tag{6}$$

where PRECIP is precipitation (mm), SURQ is surface runoff (mm), GWR is groundwater recharge (mm), AET is actual evapotranspiration (mm), SW is soil water storage (mm), and n is land use type. Water yield is also an essential WRES indicator related to the sustainable management of water resources in the study area and the key to river regime sustainability. Water yield is the amount of water from each ecosystem that enters the water body. Water yield has a complex component consisting of surface runoff with a short concentration time and lateral flow and baseflow with a longer concentration time.

$$\text{WYLD}\_{\text{n}} = \text{SURQ}\_{\text{n}} + \text{LAT}\_{\text{n}} + \text{BFO}\_{\text{n}} \tag{7}$$

where WYLD is water yield (mm), SURQ is surface runoff (mm), LAT is lateral flow (mm), and BFO is baseflow (mm). This study also assesses WRES on an annual and seasonal scale. Temporal assessment is significant to rationally allocate water resources, especially

in areas with seasonal excess water and drought. SWAT output is separated based on the monthly rainfall pattern in one year, the wet month when the rainfall is >200 mm, and the dry month when the rainfall is <100 mm.

#### **3. Results**

#### *3.1. Soil Water Retention Due to Soil Compaction*

Bulk density and porosity as indicators of soil compaction were observed in the topsoil because agricultural activities that encourage soil compaction occurred in the topsoil compared to the subsoil. The results showed that the total pore size had the following trends: FP (55.6%) > YOP (52.0%) > AGF Ex-MOP (49.3%) > MOP (49.0%), while the bulk density trends as follows: MOP (1.35 g/cm3) > AGF Ex-MOP (1.34 g/m3) > YOP (1.12 g/cm3) > FP (0.91 g/cm3). The relationship between bulk density and soil porosity is reciprocal, meaning that higher bulk density will reduce porosity and vice versa. FP has the lowest bulk density and highest porosity, while MOP has the highest bulk density and lowest porosity. Based on ANOVA and DMRT, FP bulk density was significantly the smallest (*p* ≤ 0.05), and soil porosity was significantly higher (*p* ≤ 0.05) than YOP, MOP, and AGF Ex-MOP (Table 1). YOP bulk density was significantly different (*p* ≤ 0.05) from MOP and Ex-AGF, but YOP porosity was not significantly different (*p* > 0.05) from FP and MOP. Meanwhile, although slight differences exist, AGF Ex-MOP and MOP bulk density and porosity were not significantly different (*p* > 0.05)

**Table 1.** The mean values of soil porosity and bulk density.


\* Significant at 95% level, n: not significant at 95% level; a,b,c the mean value followed by the same letter does not differ according to the DMRT.

Due to soil compaction, bulk density and soil porosity changes lead to soil water retention changes. The soil water retention curve (SWRC) is presented in Figure 2, showing that the trend of soil water retention in each land use has the same pattern as the soil porosity: FP > YOP > AGF Ex-MOP > MOP. These results indicate that FP has the highest soil water retention for the same potential energy, while MOP has the lowest. The best-fit van Genuchten parameters obtained from adjustments are presented in Table 2. The SWRC of each land use has high suitability to the observed data, evidenced by a low RMSE and high R-squared. Parameter n and m are related to the slope of the retention curve after the inflection point (S) [10]. The S value is often used to describe the level of soil degradation. The results obtained indicate that MOP has the highest level of soil degradation, while FP is the lowest.

**Table 2.** SWRC best-fit parameters and S value.


**Figure 2.** Soil water retention curve for each land use.

Precipitation or irrigation that infiltrates the soil column moves freely under gravitational force, and the soil matrix binds the rest by adhesive force. Gravitational water occupies drainage pores or the range of soil pores between saturated water content and field capacity (pF 0–pF 2.54). FP has higher saturated water content, field capacity, and drainage pores (Table 3), which indicates that FP has more gravitational water than YOP and MOP. Land use with high gravitational water has implications for higher percolation, lateral flow, and groundwater recharge, evidenced by SWAT simulation. Soil water bound by soil matrix is divided into available water content (AWC) or water that plant roots can still absorb and permanent wilting point (PWP) or water that plant roots can no longer absorb. AWC is the water content occupying the available water pore space (pF 2.54–pF 4.2). In contrast, PWP occupies the unavailable water pore space (pF ≥ 4.2). FP with high available water pore space implies that FP holds more soil water as a source for the plant uptake and evapotranspiration process.


**Table 3.** Soil water retention characteristic on the same potential for each land use.

#### *3.2. Evaluation of Water Regulation Ecosystem Services*

The performance of the calibrated SWAT was evaluated quantitatively based on statistical values compared to the criteria recommended by Abbaspour [27], Moriasi [29], and Ayele [30]. The model's performance is very good for the calibration period, with NSE 0.78 and R-squared 0.88, and suitable for the validation period, with NSE 0.67 and R-squared 0.83. In addition, the p-Factor and r-Factor obtained are 0.93 and 1.09 for the calibration period and 0.75 and 0.54 for the validation period, which indicates model uncertainty is acceptable. Based on the literature review, 20 key parameters capable of capturing the main WRES were selected for calibration and 7 of them were the most sensitive parameters based on the Latin hypercube sensitivity analysis (Table 4).


**Table 4.** SWAT calibrated parameters with their range and best-fit values.

<sup>a</sup> v: replace the initial value with the best fit value, r: multiply the initial value with the best fit value; <sup>b</sup> parameter is sensitive when *p*-value < 0.05 or |T-stat| > Tα. df, sensitive parameters is marked with (\*).

All calibrated parameters are CN2 (curve number in average moisture conditions) and OV\_N (manning "n" coefficient for overland flow) that related with surface runoff; LAT\_TTIME (lateral flow travel time) that related with lateral flow; CH\_N2 (manning "n" coefficient for the main channel), CH\_K2 (hydraulic conductivity of the main channel), CH\_N1 (manning "n" coefficient for tributary channel), and ALPHA\_BNK (riverbank recession constant) that related with streamflow routing; CANMX (maximum canopy storage), ESCO (soil evaporation compensation coefficient), and EPCO (plant uptake compensation factor) that related evapotranspiration; ALPHA\_BF (baseflow recession constant) and GWQMN (baseflow threshold) that related baseflow; and REVAPMN (water level threshold for "revap"), GW\_DELAY (groundwater delay), and RCHRG\_DP (deep aquifer recharge proportion) that related groundwater recharge. Furthermore, all observed parameters are SOL\_BD (bulk density), SOL\_Z (soil depth), SOL\_AWC (available water content), and SOL\_CBN (soil carbon content) that related with soil water storage; and SOL\_K (soil permeability) that related with lateral flow. The initial value of curve number (CN2) and manning "n" coefficient for overland flow was based on frequently used and reliable literature.

The three parameters related to streamflow routing are sensitive parameters, such as ALPHA\_BNK, CH\_K2, and CH\_N2. ALPHA\_BNK is the most sensitive parameter, indicated by the highest |T-stat|. The sensitivity of ALPHA\_BNK, CH\_K2, and CH\_N2 indicates that the flow routing mechanism influenced by these parameters greatly determines the streamflow dynamics. The sensitivity of CH\_K2 shows that streamflow is strongly influenced by two-way interactions between rivers and shallow aquifers, where this interaction only occurs in intermittent rivers. Rivers receive water from shallow aquifers when the water table level exceeds the riverbed (rainy season) and lose water when the water table level is less than the riverbed (dry season). The velocity of streamflow filling and loss is closely related to the hydraulic conductivity of the soil layer (CH\_K2) between the riverbed and the shallow aquifer. CH\_N2 is also a sensitive parameter, which means that the flow velocity greatly determines the streamflow dynamics. Higher CH\_N2 is associated with lower flow rates, while lower CH\_N2 is associated with higher flow rates.

The other sensitive parameters are CN2, ESCO, SOL\_BD, and GWQMN. As a component that dominates streamflow during the rainy period, the magnitude of surface runoff

calculated from the CN2 value also dramatically determines the streamflow dynamics. Previous studies have shown that CN2 is always a sensitive parameter when the CN method is chosen for rainfall–runoff modeling [33]. CN2 has a value range of 0 to 100, but the often-used values are in the range of 25 to 98. The greater the CN2 value, the higher the surface runoff generated from rainfall. SOL\_BD is a parameter that determines the soil water retention, which then implies soil moisture dynamics. Soil moisture dynamics are necessary to determine surface runoff, lateral flow, groundwater recharge, and actual evapotranspiration. The last, GWQMN, is a parameter that determines the amount of baseflow, where baseflow only appears if the water table exceeds the GWQMN value.

The SWAT model used for evaluating WRES has been updated by adding soil water retention characteristics in the HRU scale. The soil water retention characteristics are inputted into soil attributes and affect the WRES value for each HRU after HRU definition. Soil water retention characteristics are related to calculating curve number retention parameter and defining soil moisture ranges that affect surface runoff, evapotranspiration, and soil water storage. Precipitation that reaches soil surface is generally allocated as surface runoff (SURQ), groundwater recharge (GWR), actual evapotranspiration (AET), and soil water storage (SW) (Figure 3). In forest patches (FP), 43% of rainfall is distributed as SURQ, 26% as AET, 25% for GWR, and 6% for SW. On the other hand, on mature (MOP) and young (YOP) oil palm plantations, rainfall is allocated as SURQ by 56% and 73%, respectively, AET by 25% and 20%, and GWR 15% and 6%. The temporal WRES of each land use is presented in Figure 4. The area of agroforestry plots on former oil palm plantations (AGF ex-MOP) are very narrow (<0.01% of the watershed area) and do not significantly affect the landscape-scale water regulation. Therefore, these agroforestry plots are not included in the SWAT simulation.

**Figure 3.** Distribution of water regulation ecosystem service elements in forest patches (**a**), mature oil palm (**b**), and young oil palm (**c**). **Notes: SURQ:** surface runoff, **AET:** actual evapotranspiration, **GWR**: groundwater recharge, **STORAGE**: soil water storage.

**Figure 4.** Hydrograph of water regulation ecosystem service elements: surface runoff (**a**), groundwater recharge (**b**), actual evapotranspiration (**c**).

Forest patches have better annual and monthly water regulation ecosystem services than young and mature oil palm plantation, evidenced by lower surface runoff, higher groundwater recharge, higher actual evapotranspiration, higher soil water storage, lower water yield in wet months, and higher water yield in dry months (Figure 5). The annual water yield of forest patches is 2222 mm, while the average water yield in the wet and dry months is 220 mm and 95 mm. On the other hand, the annual water yield of young and mature oil palms were 2465 mm and 2298 mm, wet month water yields were 261 mm and 237 mm, and dry month water yields were 55 mm and 72 mm. The annual and seasonal water yield difference between forest patches and young oil palms is very high and decreases as the oil palm grows [2]. An increasing proportion of surface runoff and decreasing lateral flow and base flow to water yield in mature and young oil palm led to more significant water yield seasonal variability despite the higher annual water yield (Figure 6). It means that oil palm water yield became concentrated in the rainy season and decreased in the dry season, while forest patch water yield is less concentrated in the rainy season and more available in the dry season.

**Figure 5.** Annual (**a**) and seasonal mean of water regulation ecosystem services (**b**) and water yield (**c**) elements.

Actual evapotranspiration (AET) is also an essential component of water regulation ecosystem services because it plays a role in crop production (water use) and microclimate regulation. The annual actual evapotranspiration of forest patches was 819 mm, with a monthly average of 77 mm in the wet months and 42 mm in the dry months. Meanwhile, the annual actual evapotranspiration of young and mature oil palms was 773 mm and

616 mm, with an average of 73 mm and 59 mm in the wet months and 34 mm and 27 mm in the dry months. Potential evapotranspiration (PET) and soil moisture at available water content (AWC) are the main limiting factors for AET. PET limits AET in saturated soil, while AWC limits AET in unsaturated soil. FP with higher AWC implied higher AET than YOP and MOP. AET values were also controlled by canopy cover, quantified by the leaf area index (LAI) by direct measurement using hemispherical photos. FP with LAI of 3.1 ± 0.63 had a higher AET than MOP with LAI of 1.27 ± 0.19. FP has an AET/PET coefficient of 0.79, while YOP and MOP are 0.59 and 0.75. This value means that, from 100% of the potential energy allocated for the evapotranspiration, FP uses 79% of that energy to evaporate water, while YOP uses 59%, and MOP uses 75% of potential energy.

**Figure 6.** Distribution of water yield elements in forest patches (**a**), mature oil palm (**b**), and young oil palm (**c**).

#### **4. Discussion**

Water regulation through the soil layer is the primary process determining how much water flows above the soil surface as surface runoff, returns to the atmosphere through evapotranspiration, percolates into the aquifer, and is stored in the soil pores. Water regulation is highly dependent on soil quality, which is indicated by soil physical, chemical, and biological characteristics that vary in response to soil type, land use, and topography [34]. Water regulation through the soil layer also determines the retention and movement of dissolved nutrients, such as nitrogen, phosphorus, and other nutrients. Therefore, good soil quality needs to be maintained to support sustainable agricultural production. However, a decline in soil quality has led to soil degradation mainly caused by anthropogenic pressures exerted on the soil beyond its carrying capacity [35].

Soil degradation in the form of soil compaction is a consequence that must be accepted due to natural ecosystem changes to agroecosystems, especially intensive monoculturecultivated agriculture, such as oil palm plantation [34,36,37]. Soil compaction in the oil palm plantation is characterized by increased soil bulk density and decreased porosity. On the other hand, ecosystem restoration, such as preserving forest patches or constructing agroforestry islands inside oil palm plantations, will improve soil structure indicated by increased porosity and reduced bulk density. Soil compaction in oil palm plantations is

mainly caused by soil structure detachment by heavy equipment pressure during soil tillage and harvesting, intensive inorganic fertilization, and decreased soil organic matter and vegetation cover [38].

Soil compactions have negative effects on the soil hydrological process, such as reduction in infiltration capacity, permeability, and soil water retention that are related with reduction in soil porosity [39]. This phenomenon, in turn, has implications for the WRES degradation. For example, [5,6] reported that the study area experienced increasing problems related to water resources due to soil compaction in the oil palm plantations area, especially the prolonged water shortage during the dry season. Based on soil water retention curve interpretation in MOP and YOP, the soil is degraded and less structured, has less organic matter, and has fewer pores. Meanwhile, FP has more structured soil with more organic matter and pores. According to [10], the SWRC slope (S value) correlates with soil bulk density, porosity, and soil organic matter content. The S value will decrease along with the increase in bulk density, decrease in porosity, and decrease in soil organic matter.

The higher the level of soil compaction causes the macropores to be further reduced, which causes the field capacity, saturated water content, and gravitational water to decrease. The decrease in gravitational water can be seen from the reduction in lateral flow, percolation, and groundwater recharge. On the other hand, soil compaction effect on water retention characteristics is almost nonexistent at very high matrix suction (pF > 4.2) [40,41]. Water retention at a high matrix suction tends to be influenced by textural pores associated with clay content. The higher the clay content, the higher the water retained by the textural pores [42,43]. Soil compaction only changes the structural pores and does not change the textural pores. For the same soil type, although there are differences in AWC and gravitational water associated with different land management, the water content in the high suction matrix tends to be the same [41]. The modeling result proves that each SWRC in the study area tends to coincide when the suction matrix gets bigger, considering that the soil in the study area has relatively the same clay content based on soil texture data from ground survey and laboratory analysis.

The SWAT model assisted in WRES upscaling from plot scale to landscape and timeseries scale. Nevertheless, before further interpretation, it is necessary to evaluate the reliability and uncertainty of the SWAT model through the NSE, R-squared, p-Factor, and r-Factor values. The objective function used during the calibration process is NSE, meaning that the value of each parameter will be optimized from its initial value until it reaches the desired NSE value. When the system calculates the NSE value, other statistical values, such as R-squared, will be adjusted automatically. In addition, evaluation of model uncertainty is vital because SWAT provides various combinations of parameters and different modules that produce one of the same outputs but produce another very different output. A good model has a p-Factor value close to 1 and an r-Factor close to 0. However, to achieve a p-Factor close to 1, it is necessary to sacrifice the r-Factor value away from 0 and vice versa. Therefore, the best simulation is defined as a balanced p-Factor and r-Factor when it reaches the highest NSE or R-squared value during the calibration and validation periods. To obtain a balanced p-Factor and r-Factor, we must arrange each model parameter's lower and upper bound.

The calibrated results are still in the reliable category, and the model's uncertainty is still acceptable, although the model's performance for the validation period is not as good as the calibration period. The high values of R2 and NSE in the calibration and validation periods indicate that the combination of calibrated parameters can capture the impact of daily meteorological input variations on daily WRES. Moriasi [29] stated that the model reliability criteria could be used for model evaluation on a monthly and daily scale. However, with the same criteria, model evaluation for the daily scale is generally more robust than the monthly one because daily output captures variations and inaccuracies arising from parameter uncertainty and daily input data (SWAT input must be daily, while the output can be daily or monthly). More accurate modeling of WRES for micro-watersheds can serve as complementary information for environmental restoration planning at a local

scale. A calibrated and validated SWAT model with acceptable reliability and uncertainty was then used to evaluate WRES of forest patches among oil palm plantations.

Due to differences in soil water characteristics, WRES varies significantly between land uses and management, so changes in land use and management will impact the change in WRES at the landscape scale by routing mechanism. Meanwhile, oil palm is the dominant vegetation in the study area, so the oil palm WRES dominates the landscapescale water balance. In addition, given the dynamic nature of oil palm plantations that are cleared and replanted every 20–30 years, it is necessary to understand oil palm WRES when they are young and mature (yielding plants). The transition from YOP to MOP occurs 8–9 years after planting, when the canopy cover reaches its maximum value. The simulation consistently shows that preserving forest patches among oil palm plantations has implications for decreasing SURQ, increasing SW, and increasing GWR. These advantages have consequences for the FP water yield, which is higher in the dry month and lower in the wet month so that water is still available in the dry season and does not overflow in the rainy season. The lower WYLD FP compared to oil palm plantations is supported by [44], which states that the response to WYLD in agroecosystems, especially oil palm plantations, tends to be higher than forest.

SURQ occurs when rainfall exceeds infiltration capacity, where lower infiltration capacity in oil palm due to soil compaction causes more rainfall to turn into SURQ. Two factors caused the higher SURQ, and more insufficient water storage in oil palm plantations: (1) decreased canopy and ground cover in YOP and (2) intensive soil compaction due to MOP harvesting. The decrease in canopy and ground cover causes the SURQ rate to be faster, time concentration to be shorter, and decreases in interception and evapotranspiration. The rough surface of the FP due to more complex canopy stratification, cover crops, and forest litter, besides lowering SURQ, also slows the SURQ rate on its way to water bodies. High infiltration capacity also increases LAT, BFO, and GWR. LAT and BFO are relatively stable WYLD parts because they have a slower rate to reach water bodies. LAT and BFO indicate the river regime's sustainability as it maintains water availability on a day without rain. An increase in SURQ and a decrease in LAT, BFO, and GWR lead to an increased risk of flooding, drought, and water scarcity. Changes in local water resources had become a significant concern for residents in the study area, including the faster shallow aquifer depletion during the dry season and high fluctuations in streamflow between the rainy and dry season [5].

High canopy cover in FP also causes soil water storage to be more compensated for the AET besides reducing SURQ. High AET increases the initial abstraction and decreases rainfall proportion for SURQ. Three factors limit AET: PET as an energy source, AWC as a water source, and stomatal conductance (correlated with LAI) as the water pump from the soil to the atmosphere. A high AET causes soil moisture to decrease faster so that the soil water changes more quickly from saturated to unsaturated conditions. Unsaturated soils have higher infiltration rates and lower SURQ rates than saturated soils, so, in this case, AET plays a role in reducing runoff. AET is also an essential element of surface energy balance related to microclimate regulation. A higher AET for the same net radiation will reduce the proportion of energy for atmosphere heating (sensible heat). The increase in air temperature in the study area is evidence of the relationship between lower AET in oil palm plantations (especially YOP) and atmospheric warming. Beside depletion in groundwater, the residents in the study area also feel the air has become much warmer since oil palm plantations have dominated the landscape [5]. In addition, the comparison of lower PET and AET in oil palm plantations (especially in the dry season) indicates that these land uses have a higher water deficit than FP.

FP has properties like a sponge and a pump in the hydrological cycle. Meanwhile, MOP only has pump properties, and YOP does not have these two properties. FP act like sponges by increasing soil water retention and releasing it slowly through LAT and BFO because it has more soil pores. Water is stored and maintained on days without rain and even remains available until the dry season through these properties. On the other hand, FP also act like pumps by evaporating a large amount of soil water into the atmosphere. The proof that MOP only has pump properties is that MOP evaporates a large amount of water, although not as much as FP, but cannot store large amounts of soil water. YOP does not have sponge and pump properties, indicated by low AET and low soil water storage. The unavoidable soil degradation due to soil compaction is the leading cause of oil palm plantations losing their sponge properties. The loss of sponge properties due to soil compaction has three consequences: a decrease in infiltration leading to an increase in SURQ, a reduction in the GWR, and a decrease in AWC leading to a reduction in AET.

The anecdotal information that oil palm plantations are "water-greedy crops" [5,45] because they have higher AET and are associated with water scarcity is inaccurate. The MOP AET tends to be less than equal to FP [2], while YOP AET is much lower. More scientific evidence for this water scarcity case is that oil palm plantations encounter soil compaction so that more rainfall flows become SURQ than stored in soil pores. Oil palm water use is greater than the stored water because of improper management, where high AET is not accompanied by high AWC and GWR, as is the case with FP. Higher AET also reduces shallow aquifers through the capillary water movement to plant roots when soil moisture is insufficient to compensate for AET, which causes the groundwater to be dwindled and become unavailable during the dry season.

WRES sustainability is mainly determined by land use and management, which changes the soil and surface and affects rainfall distribution into SURQ, GWR, and AET [46]. According to locals' information that there is a groundwater decline during the dry season as the primary water source, the desired alternative for water management is to reduce SURQ and increase GWR. Landscape SURQ can be reduced and GWR can be increased through a multifunctional landscape, i.e., retaining the remaining FP among oil palm plantations. By maintaining FP or agroforestry as a high conservation area with an optimal area, soil and water conservation becomes a more suitable alternative for improving landscape WRES while maintaining oil palm productivity. FP with good soil structure and more complex canopy stratification enhance the WRES, so their existence in a landscape dominated by oil palm plantations is vital to maintaining the watershed's ecological integrity. In addition, the presence of FP inside oil palm plantations can synergize the provisioning and regulating services, where oil palm provides provision ecosystem services and FP provide regulation ecosystem services. The recommended multifunctional landscape is to maintain FP or create agroforestry islands separately around oil palm plantations. If oil palm is mixed with forest vegetation, there will be competition for light and water, which hinders the productivity of the entire vegetation.

Multifunctional landscapes also enhance other ecosystem services besides WRES, such as biodiversity conservation and erosion prevention [47]. Multifunctional landscapes are also an effort to adapt to the negative impacts of climate change. Changes in rainfall patterns due to climate change, where rainfall is predicted to increase in wet months and decrease in dry months causes WRES variations to become more extreme [48]. Further research is needed regarding the spatial configuration of multifunctional landscapes with the most optimum ecosystem functions and economic benefits based on ecosystem services trade-offs. Another aspect that needs to be considered is that the chosen multifunctional landscapes policies must be site-specific and examine the local biophysical characteristics of the landscape, such as soil type, topography, and elevation [49,50]. Even more broadly, they need to include social and economic factors. For example, soil type in the study area belongs to the hydrologic soil group (HSG) C based on the observation of soil permeability, where the runoff coefficient of each land use is higher than HSG A or HSG B. In addition, the topography in the study area belongs to the relatively flat areas, where the runoff coefficient is lower than the steeper slope. Environmental planners require these site-specific quantitative relationships to balance landscape-scale ecological and socio-economic functions.

#### **5. Conclusions**

Multifunctionality landscapes approach through maintaining forest patches between oil palm plantations, can improve landscape WRES, shown by a decrease in surface runoff, an increase in groundwater recharge, an increase in soil water storage, and an increase in actual evapotranspiration. As a result, water is not concentrated in the rainy season and remains available in the dry season. The forest patches can improve landscape WRES because they have good soil hydrological characteristics, so the existence of forest patches is essential to maintain the ecological integrity of the watershed. Soil hydrological characteristics in forest patches are indicated by the lowest bulk density and the highest soil porosity compared to oil palm plantations. Good soil hydrological characteristics have implications for increasing soil water retention, imply more soil water is stored and available to plants (available water content), and flows through soil pore spaces to fill aquifers and river networks (gravitational water). Meanwhile, soil compaction increases bulk density, decreases porosity, and decreases soil water retention in oil palm plantation.

The calibrated SWAT is reliable and acceptable model, shown by the high NSE and R-squared value and balanced p-Factor and r-Factor values. The SWAT model consistently proves that forest patches have sponge and pump properties in the hydrological cycle. The sponge properties are related to the optimal distribution of soil porosity so that water is stored in the rainy season and flows slowly during the dry season. Meanwhile, the pump properties are related to plant tissue, which plays a role in absorbing and returning soil water to the atmosphere to maintain microclimate stability. Mature oil palm can only evaporate large amounts of water like forest patches but cannot retain some soil water because of degraded soil pores due to soil compaction. The pump properties, which are not accompanied by the sponge properties, cause the water use to be greater than the stored water. Therefore, the multifunctional landscape approach by conserving forest patches between oil palm plantations is one approach that can improve the sustainability of oil palm plantations. The multifunctional landscape can synergize provisioning services from oil palm plantations and regulating services from forest patches.

**Author Contributions:** Y.K. and S.T. designed the research. Y.K. wrote the manuscript. S.T., T.J., E.D.W. and B.S. reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by PMDSU (Program Pendidikan Magister Menuju Doktor untuk Sarjana Unggul) scholarship from Ministry of Education, Culture, Research, and Technology, Republic of Indonesia.

**Data Availability Statement:** The datasets presented in the study are included in the article material; further inquiries can be directed to the corresponding author/s.

**Acknowledgments:** This research was supported by CRC-990 EFForTS (Collaborative Research Center 990 Ecological and Socio-economic Functions of Tropical Lowland Rainforest Transformation System) in form of access to the study site and provision of daily meteorological data.

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

#### **References**


**Aikaterini Voudouri, Evgenia Chaideftou and Athanassios Sfougaris \***

Crop Production and Rural Environment, Laboratory of Ecosystems and Biodiversity Management, Department of Agriculture, University of Thessaly, Fytokou Str., N. Ionia, GR-38446 Volos, Greece; k.voudouri@hotmail.com (A.V.); eugeniachd@gmail.com (E.C.)

**\*** Correspondence: asfoug@agr.uth.gr

**Abstract:** The topsoil seed bank was studied in four types of agricultural bird habitats: fields with cereals, maize, clover and tilled fields of a Mediterranean plain to determine the potentially richest habitat based on food supply for the wintering farmland birds. The diversity and abundance of topsoil seeds differed between seasons but did not differ significantly between habitats. The cereal habitat was the richest in food supply for the overwintering of farmland birds. The topsoil seed bank was dominated by *Chenopodium album*, *Polygonum aviculare* and *Amaranthus retroflexus*. The findings of this study provide insight for low-intensity management of higher-elevation mount agricultural areas of southern Mediterranean by preserving seed-rich habitats for farmland avifauna.

**Keywords:** topsoil seed bank; farmland bird diet; agricultural ecosystem; biodiversity; habitat

#### **1. Introduction**

The management of agricultural land has greatly changed over recent decades. This has resulted in different physiognomy and a reduction of agricultural biodiversity and heterogeneity [1–4]. The depletion of the natural transient soil seed bank during cultivation was one of the changes (e.g., [5]). Shifts in agricultural management also led to decline of rural birds [6–13].

Mosaics of low-intensity cultivation in the Mediterranean areas may preserve high diversity of bird species, but intensification or land abandonment probably do not benefit biodiversity [4,14]. Reviews identified that agricultural intensification [15] and concomitant abandonment [16] remain the major threat to agricultural ecosystems of the 21st century across Europe and elsewhere (e.g., [17]) with many ecological and biodiversity impairments. As a consequence, investigation of floristic and seed diversity and abundance, along with the physiognomy of the rural landscape, is necessary for identifying the most interesting bird habitats. These habitat features, i.e., high quality or food resources or aboveground floristic components like stubbles or semi-natural with natural habitats suitable for breeding (low intensity farmland with steppe-like vegetation), which can be proved to be beneficial for birdlife per case facilitate biodiversity maintenance [6,18].

Approximately 30% of the bird species being "Species of European Conservation Concern" exploit agricultural ecosystems [19,20]. In rural landscape across Europe, food availability (plant and seed food items) is reduced especially during winter [3,21], and nesting habitats in spring are deficient for many bird species that have declined over recent decades [22]. For instance, seed-eating birds face the risk of limited accessibility to preferred seeds when vegetation is dense in uncultivated areas close to farmlands [23] or they feed in stubble of reduced quality due to modernized techniques in cereals harvested as arable silage [24].

Especially in winter the most important food resource in the rural landscape for birds is the soil seed bank (e.g., [3,25]). However, research has been focused on the impacts of agricultural practices on the seed bank composition, abundance and vertical distribution of

**Citation:** Voudouri, A.; Chaideftou, E.; Sfougaris, A. Topsoil Seed Bank as Feeding Ground for Farmland Birds: A Comparative Assessment in Agricultural Habitats. *Land* **2021**, *10*, 967. https://doi.org/10.3390/ land10090967

Academic Editors: Chiara Piccini and Rosa Francaviglia

Received: 13 August 2021 Accepted: 3 September 2021 Published: 14 September 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/).

weeds in the rural landscape regarding their persistence (e.g., [26]), fertilization (e.g., [27]) crop rotation or varied tillage systems (e.g., [28,29]). As a result, the seed bank fraction on topsoil of agricultural habitats that serves as food resource to bird seed-eaters is rarely, if at all, studied.

Food resource provision [30] is crucial for wintering seed-eaters especially them of a high conservation interest. However, it is less studied in the Mediterranean regions compared to northern Europe. For conservation of bird populations which are exclusively or partially seed foragers, such as *Passer montanus* (Eurasian tree sparrow), *P. domesticus* (house sparrow), *Fringilla coelebs* (common chaffinch), *Carduelis chloris* (European greenfinch), *C. carduelis* (European goldfinch), *Miliaria calandra* (corn bunting), *Turdus merula* (common blackbird) and *Emberiza* spp. (bunting birds), a precise "instruction" of the most proper habitat, i.e., seed-rich in winter [31], is not yet defined [32].

A seed-eating bird may have preference on the seeds of particular plant species for their diet [3,22] thus seed consumption can cause seed limitation of these particular species. As a result, determination of the value a habitat carries to supply winter food to birds [33] is critical in decision-making for maintaining diversity in agricultural ecosystems [21,34,35]. Therefore, there is scope for further consideration of how we manage areas of former traditional low-intensity agriculture in Europe, supported by European subsidies [24], and principally seed-rich habitats such as cereal stubbles in certain seasons of the year (e.g., [36,37]).

This study aimed at determining the potentially richest habitat (for food supply) for the wintering of rural avifauna in the Dolichi plain of Elassona region. To this aim, the study investigated the effect of the type of crop (habitat for birds) on the topsoil seed bank in four arable fields undergoing post-dispersal consumption of seeds by farmland birds. This topsoil seed bank entails a fraction that serves as food source especially to the seed-eaters.

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

#### *2.1. Study Area*

The study area covers an area of 40 km2, where the settlement of Dolichi and the municipality of Livadi are situated. Livadi is located at an altitude of 1100 m. Dolichi lies 5 km from the foothills of Mount Olympus and 21 km from Elassona, at an altitude of 590 m asl. The inhabitants of Dolichi, numbering 473 (in 2001) are principally involved in land cultivation and animal husbandry.

The settlement of Dolichi is located in the center of a cultivated plain and is surrounded by hilly and mountainous natural ecosystems. The main crops are cereals, with wheat (*Triticum* sp.) dominating over the barley (*Hordeum* sp.). Previously, the second largest crop was tobacco, but due to the regime of European subsidies tobacco cultivation was abandoned, therefore the second rank is currently held by maize crops (*Zea mays*). Legumes and vegetables are cultivated on a smaller scale. Clover (*Trifolium* sp.) in particular is cultivated in a mixed agricultural-livestock farming system. In 2005, vetch crops exceptionally dominated due to compulsory crops rotation in line with the Codes of Good Agricultural Practice. Ecosystems with a high plant cover such as grasslands, hedgerows, uncultivated vegetation strips and riparian zones are largely present in the area. The cultivation system is mechanized, but clearly less intensive than the one of the main plain of Thessaly. In the area many plantations of locust (*Robinia pseudoacacia*) have been established through subsidies of the agri-environmental measure "Afforestation of agricultural land".

The landscape characteristics in this area were appropriate for the research. Approximately two thirds (2/3) of the study area are covered by anthropogenic ecosystems and only one third by natural ecosystems. Among natural ecosystems, grasslands (including ecotones) hold the largest area. Cereals comprise 90% of anthropogenic ecosystems while legumes only 5%. The remaining 5% is covered by other types of agricultural ecosystems.

In this area of study, fields are undergoing post-dispersal consumption of seeds by farmland birds, having thus a potential as feeding sources to them, with non-cultivated and cultivated plant species ([38]; Table A1 of Appendix A). The study took place in four selected habitat types which actually are arable cropped fields: cereals, maize, clover and tillage, with average surface area of 2.35 ± 0.3 ha ([38]; Table A2 of Appendix A). According to [38], within this area the above-ground non-cultivated species richness and (%) plant cover differed from the respective ones of plant species serving as food sources to farmland birds among the four studied habitats (Table A2). This is not true though for the field physiognomy (see Table A2).

The avian diversity in this area was also proper for the research aims. Overall, 33 bird species were recorded in the study area (unpublished data, see Table A3 in Appendix A). A high majority, that is, 26 bird species comprising 79% of the recorded bird species in the study area, are classified as seed-eaters and all birds listed in Table A3 are present in all studied habitat types (namely the crop types: cereals, maize, clover and tillage).

#### *2.2. Research Design*

The sampling area is shown in Figure 1. Four types of fields were sampled and analyzed between September 2006 and March 2007 in the current study: cereals, maize, clover and tillage (bare soil during winter). Since these fields include 'micro-sites' providing refuge and food resources to rural bird populations, they are referred to here as "habitats". There were crops and stubbles in the maize and cereal fields during winter.

**Figure 1.** Map of Dolichi plain in the area of Elassona (study area). Image and photo sources: The blank map of Greece (on the top left) is by Lencer, CC BY-SA 3.0, https://commons.wikimedia.org/ w/index.php?curid=4432468. The satellite map image is extracted and edited from Google Earth (https://earth.google.com/) on 22 November 2017. The landscape photo of the Dolichi plain is taken from [38].

Plots of approximately 20 m2 each were randomly chosen so that the major species are represented in the cultivated area. A total of 36 plots, in total 846.7 m2, were recorded and sampled. The number of studied plots (n) per crop type was: n = 10 in winter cereal fields (197.5 m2), n = 10 in maize fields (274.5 m2), n = 6 in clover fields (171.5 m2) and n = 10 in tillage (203.2 m2) (Table A4).

#### *2.3. Topsoil Seed Bank Sampling, Seed Extraction and Identification*

During fall-end and winter-start of 2006, within 21 plots (cereal n = 6, maize n = 6, clover n = 3 and tillage n = 6) randomly selected out of the total of 36, soil cores were sampled (soil corer of 15 cm diameter, and 1 cm depth). In randomly selected plots (at every second measurement of plant cover), 10 soil samples were collected (R = 10) across the

diagonal of the plot and were placed into encoded polyethylene bags that were transferred to the laboratory. The second soil core sampling took place in spring of 2007 following exactly the same protocol, though the number of plots differed due to the seasonal change of landscape in Dolichi plain. A total of 12 plots were totally studied in spring: cereal (n = 4), maize (n = 2), clover (n = 3) and tillage (n = 3). Seeds up to one centimeter of soil depth have been sampled. Therefore, the seeds and seed bank are referred in this article as 'topsoil'.

Soil core samples were retained at 4 ◦C for 24 h to avoid seed germination. Then seeds were isolated from soil phase using sieves and were identified at species level in petri dishes under stereoscope and magnifying lens using a series of keys (Appendix B). In addition, plant specimens were collected in spring and autumn of 2006 so that all plant species the seeds of which are potentially present in the topsoil seed bank of autumn 2007, are included in a reference plant collection that facilitated the seed identification (sources are listed in Appendix B). Classification of plant species on the basis of their significance to the farmland bird diet is presented in Tables A1 and A5 of Appendix A.

#### *2.4. Data Analysis*

The species richness of the seed bank was estimated as the number of species per m2. The Shannon index (entropy) was also estimated for the topsoil seed bank. The seed abundance of the topsoil seed bank was estimated as the average number of seeds per square meter (m2).

The Shannon index (entropy) and the seed abundance of the topsoil seed bank were tested by estimating the differences between habitats (cereal, maize, clover and tillage) and seasons (winter and spring) using generalized linear mixed (GLM) effect models. For Shannon index, the GLM for Gaussian family with random intercept of plot were used. For seed abundance, the GLM for Poisson family (with log link function) with random intercept of plot were employed. For model selection, model with and without interaction between season and habitat were compared with simple generalized linear model (for Gaussian family in the case of Shannon Index; for Poisson family (with log link function) in the case of seed abundance). In each case, all three models were compared using *p*-value of ANOVA and Akaike Information Criterion (AIC). Signs of heteroscedasticity (residuals vs. fitted plot) and normality of residuals (q-q plot) were also tested for identifying the best performing model. Overdispersion of the Poisson model was also checked and, if needed, the analysis was reconducted using generalized linear (simple or mixed effect) model with negative-binomial family with log link function.

Predictions were generated with and without inclusion of random effects. The 95% confidence intervals were estimated with bootstrapped simulation (n = 1000) using the bootMer function. Post-hoc (Tukey all-pairs) comparisons were conducted.

The data were processed in R 4.1.0 [39] using the packages: broom.mixed [40], dplyr [41], ggplot2 [42], lme4 [43], lmerTest [44], MASS [45], multcomp [46], and tidyr [47].

#### **3. Results**

#### *3.1. Composition of the Topsoil Seed Bank as Food Source to Farmland Birds*

Overall, the soil seed bank sampling and seed identification resulted in a total of 66 plant species, 64 of which are non-cultivated and the other two are the cultivated species *Triticum aestivum* and *Zea mays*. Out of these 66 plant species, 49 species were classified as having a level of significance as food items to the farmland birds (Table A5).

Soil seed bank sampling during winter in 36 fields resulted in a total of 62 plant species, 15 of which were present in all habitats while the soil seed bank sampling during spring in 12 fields resulted in 39 identified species. A total of 26 and 21 species of seeds serving as food source to farmland avifauna were identified in winter and spring topsoil seed bank, respectively (Table A5).

The highest number of species with significance as food sources to birds was recorded in cereals (34 species) while the lowest in tillage (23 species) (see Table A5). The habitats can be ranked on the basis of the number of significant species in serving as food sources to birds as follows: tillage (23) < clover (25) < maize (31) < cereal (34).

The winter seed bank was dominated in all studied habitats by *Chenopodium album*, *Polygonum aviculare* and *Amaranthus retroflexus*. The last two species also dominated the spring seed bank. Apart from these two species, the following species predominated in the spring seed bank: *Lithospermum arvense* in cereals, *Amaranthus albus* in maize, *Echinochloa crus-galli* in clover and *Digitaria sanguinalis* in tillage (Table A5).

Commonly, 35 species were present in both spring and winter while only seven species (*Brassica juncea*, *Sinapis arvensis*, *Silene dioica*, *Chenopodium album*, *Polygonum aviculare*, *Portulaca oleracea*, *Amaranthus retroflexus*) were present in all habitats in both seasons (Table A5).

These common and dominant plant species, with the exception of *Amaranthus* sp. and *P. oleracea*, are classified to highest significance as food sources to farmland birds (Table A5).

#### *3.2. Shannon Entropy and Seed Abundance of the Topsoil Seed Bank*

#### 3.2.1. Model Selection

Mixed effect model without interaction did not differ significantly from model with interaction for both seed abundance and Shannon index (entropy). However, it differed in both cases significantly from model without random intercept. The model without interaction also performed best in terms of AIC in the case of seed abundance, and similar to model with interactions in the case of Shannon entropy. Visual inspection showed no clear signs of heteroscedasticity nor deviation from normality of residuals.

However, in the case of seed abundance, the model showed significant overdispersion (dispersion ratio = 12.97, Pearson's Chi-Squared = 2931.16, *p* < 0.001). Thus, instead of using Poisson family model, negative binomial family was used for model estimation. The new model also performed better in case of AIC (Tables 1 and 2). Visual inspection showed that normality of residuals was slightly worse than in Poisson family model.

**Table 1.** Generalized Linear Model comparison for Shannon entropy and seed abundance in the topsoil seed bank.


+ indicates it is additive model (y = spring + winter + tillage + clover + maize + cereals). \* indicates there is interaction between habitat and season (y = spring + winter + tillage + clover + maize + cereals + spring \* tillage + spring \* clover + spring \* maize + spring \* cereals + winter \* tillage). \*\* indicates statistical significance, *p* < 0.05.


**Table 2.** Generalized Linear Model comparison for Shannon entropy and seed abundance in the topsoil seed bank on the basis of AIC.

+ indicates it is additive model. \* indicates there is interaction between habitat and season.

The calculated parameters of the model for Shannon entropy and seed abundance are summarized in Table 3.

**Table 3.** Model summary for Shannon entropy and seed abundance.


3.2.2. Shannon Entropy of the Topsoil Seed Bank

The Shannon entropy of the topsoil seed bank did not differ significantly between seasons or between habitats (Table 4).

**Table 4.** Tukey all-pairs comparisons for Shannon entropy of the topsoil seed bank.


In spring, the estimated Shannon entropy was higher in cereal than in the other three habitats; while in maize it was the lowest (Table 5). However, these differences were not statistically different.


**Table 5.** Model predictions compared to Shannon entropy (mean and median) of the topsoil seed bank.

Table 5 presents the model predictions, both unadjusted (not including random effect) and adjusted (including random effect), compared to the original-data mean and median.

3.2.3. Seed Abundance of the Topsoil Seed Bank

Significant differences in the seed abundance of the topsoil seed bank were observed only between seasons (winter-spring; see Post-Hoc comparisons in Table 6). In tillage the seed abundance was lower than in the other three habitats, and this difference was more pronounced with clover (Table 7). It is noted, however, that these differences in the seed abundance between habitats were insignificant and this is also confirmed by the Post-Hoc comparisons (Table 6).

**Table 6.** Tukey all-pairs comparisons for seed abundance of the topsoil seed bank.


\* indicates statistical significance, *p* < 0.05.


Table 7 shows the model predictions, both unadjusted (not including random effect) and adjusted (including random effect), compared to the original-data mean and median.

#### **4. Discussion**

#### *4.1. Effect of Habitat (Crop) Type on the Topsoil Seed Bank*

The differences in Shannon entropy of the topsoil seed bank were insignificant for the tested samples of this study. The effect of habitat (crop) type in such landscapes require, to our view, further coordinated research that would also include samples of different size, and consideration of a soil-property matrix [29].

The winter and spring Shannon entropy were lower in tillage though not significantly. This difference, although insignificant, could be explained by the widespread practice of autumn tillage which buries the surface seeds [48] thus reducing their availability in fields [49]. Tillage techniques prevent vegetation growth, seed germination and seedling growth [50] and temporarily enrich the topsoil with seeds [51]. Topsoil seeds are easily depleted from the soil surface also because they are consumed by high numbers of birds using stubbles as feeding grounds [52]. Moreover, in all studied fields where tillage was employed the seed abundance of dominant species was very high, as tillage decrease seed diversity [53]. [29] stressed the significance of such practices to preserve biodiversity in crop fields, and the complexity of it, as continuous tillage was found to have increased the soil seedbank diversity and density under specific soil conditions.

The winter topsoil seed bank is dominated only by *Chenopodium album*, *Polygonum arviculare* and *Amaranthus retroflexus* that have long-lived seeds according to [54]. By contrast, the spring topsoil seed bank reveals a different picture since apart from *Chenopodium album*, and *Amaranthus retroflexus*, four other species are prominently present: *Lithospermum arvense*, *Amaranthus albus*, *Echinochloa crus-galli* and *Digitaria sanguinalis* in cereal, maize, clover and tillage respectively. The above species, apart from *Lithospermum arvense*, are also reported to have seeds of high longevity [54]. However, note that the seeds of *Echinochloa crus-galli* found in clover and *Digitaria sanguinalis* found in tillage in our case, are classified as very-short lived for untilled systems by [55].

Since in this study only the topsoil seed bank has been investigated, no conclusion on seed persistence per species under heavy disturbance can be given. Regarding seed availability as food sources to farmland birds it should be considered that in no- or lowtillage fields where the soil is not heavily disturbed, seed predation is enhanced [56,57].

In this study, the genus Amaranthus dominated [58] maize crops. *Chenopodium album* dominated the topsoil seed bank of all habitats either in spring or in winter. This is consistent with the findings of [59], who detected high seed abundances of this annual broadleaved species in the upper 5 cm of soil irrespective of barley tillage treatments in Alaska, as well as of other authors [58,60]. *Polygonum arviculare* was dominant in winter in the topsoil seed bank, implying that autumn tillage did not bury the bulk of its surface seeds. [59] detected higher seed density of *Polygonum arviculare*, only for medium-intensity tillage treatments (disc once) in the upper 5 cm of soil.

Ref. [61] found in fields of Poland that the base of the winter diet of reed bunting *Emberiza schoeniclus* are seeds of *Chenopodium album*, *Amaranthus retroflexus*, *Setaria viridis, Stellaria media* and *Fumaria officinalis*. It should be underlined that the aforementioned differences in dominant seeds are consistent with the high spatial variability of seeds predated, such as *Chenopodium album*, given that some birds of the study area may count on alternate food resources, have preference to specific species and respond differently to different plant cover [62] and landscape composition in winter [35].

Species present in the aboveground flora and linked to disturbance in agricultural soils are *Amaranthus* spp., *Chenopodium album* and *Echinochloa crus-galli*, while *Digitaria sanguinalis*, of which seeds were dominant in clover, are linked to undisturbed soils [54] and these species have over 3-year seed longevity [54].

#### *4.2. Agricultural Habitats with a Topsoil Seed Bank Serving as Food Source to Farmland Birds*

A total of 26 and 21 species of seeds serving as food source to farmland avifauna were identified in winter and spring, respectively. Differences between spring and winter seed abundance are mostly attributed to seed consumption of species with high significance to farmland bird diet in this study. The farmland birds use stubbles more frequently in winter [63,64]. The main food source of seed-eating birds during winter is the soil seed bank [3,22,25,65].

Cereal seeds would rather show a higher potential to positively impact rural bird diversity in the studied landscape, while the structural characteristics in clover habitats might also favor birds' presence, but these are objectives of future study in a more systems' thinking approach, beyond single-farm scales [21]. [66] underlined the importance of features like hedgerows in diversifying habitats associated with many farmland bird benefits. [67] proposed that the best option for birds in winter are the seed-rich habitats while in the summer structurally and floristically rich habitats. The results of this study would rather support the findings of [18] that highlighted the importance of the presence of suitable breeding habitats in mixed landscape for farmland birds. Furthermore, in our case, whether differences between spring and winter seed abundance are attributed to seed consumption of species with high significance to farmland bird diet needs further investigation. From this viewpoint, more thorough investigation of the relationship between the richness and abundance of bird fauna and the respective parameters of seeds is necessary in the future.

Seed-eating birds are important topsoil-seed consumers inferring quantitative and qualitative changes in soil seed bank, especially in winter when plants serving as food sources to avifauna are highly reduced, compared to spring in the same fields [65]. Conversely, seed predators also have a determining role in plant communities at landscape level by impacting the abundance of specific plants of their preference, thus inferring floristic variations even at areas that are distanced in the landscape [68]. Consequently, it could not be disregarded that reduced seed availability can be a limiting factor to wintering birds, a fact that highlights the importance of interspecific competition of avian communities [69]. As such, neither the preference of seed foragers for seeds of varied seed sizes of specific annual plant species at landscape patches is to be overlooked in current and future agro-ecological management [70] nor the importance of the minimum distance of available food-resource patches in the rural landscape [31].

In this respect, as [71] supported, the landscape heterogeneity may benefit generalist birds but may mean habitat loss and fragmentation for specialists, and therefore management should not include unique standalone measures. Fragmentation and land-use changes in rural landscape also influence the soil seed bank in terms of size and composition [72]. These, and the current study findings, highlight the importance of habitat provision for farmland birds during winter and breeding seasons [21].

**Author Contributions:** Conceptualization, A.V. and A.S.; methodology A.V. and A.S.; software, A.V. and E.C.; validation, A.V., E.C. and A.S.; formal analysis, E.C.; data curation, E.C.; writing—original draft preparation, A.V., E.C. and A.S.; writing—review and editing, E.C. and A.S.; visualization, E.C.; supervision, A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Data available on request from the authors.

**Acknowledgments:** The authors would like to particularly acknowledge the contribution of Martin Hermy, at KU Leuven, for his critical review comments on a previous version of this manuscript. Special appreciation to Michał Czyz, at Coding Manatee Ninja, for advice and verification of the ˙ data analyses conducted. We also wish to thank the reviewers for their improvements on a previous version of this manuscript.

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

#### **Appendix A**

**Table A1.** Presence (+), total species richness and phenology of 61 herb-layer non-cultivated species of 25 families (11 grasslike species and 50 broadleaf (4 legumes and 46 forbs) herb-layer species) in Cereals, Maize, Clover, and Tillage in Dolichi plain during the growing season of 2006 (from [38]).



**Table A1.** *Cont.*

\* Phenological classes of herb-layer species: according to life cycle: A = Annual, B = Biennial, P = Perennial; according to biological group (BG): G = Grass, L = Legume, F = Forb, B = Bulbous (geophyte). Classification of plant life form in line with [73] in parenthesis: The = Therophyte; Hem = Hemicryptophyte; Pha = Phanerophyte; Cha = Chamaephyte; Cry = Cryptophyte; Cli = Climber. G: Cyperus is considered grass-like. Related background references: [54,74–77]. In **bold** are shown species of (least to highest) significance as food items to rural birds, while significant and highly significant species are besides underlined; classification followed [22,78,79], as well as field observations.

> **Table A2.** Mean (±Standard Deviation) above-ground variable estimations in the studied habitats [38]. Plant cover was recorded following [64]. The Field Physiognomy Index estimation followed [80].


<sup>1</sup> Species richness of non-cultivated species and of species serving as food items to birds differ with habitat (physiognomy), with the highest values in cereals (1-way ANOVA; LSD; *p* < 0.001; unpublished data from [38]).

Recordings were conducted using quadrat (1 m2) at sampling points of a plot diagonal (see also design in Table A4). The recorded variables were: 1. The total number of noncultivated plant species. 2. The percentage of the sampling plot area (1 m2) that is covered by non-cultivated plant species (%) percent of non-cultivated plant cover). 3. The number of plant species serving as food items (namely species producing seeds where birds feed on) to rural bird (food resources); this classification was based on [22,78,79] and observations in the field. 4. The (%) percentage of food items from the total surface non-cultivated plant cover serving as food items to rural birds (%) percent of non-cultivated plant cover serving as foot item to birds. 5. The (%) percentage covered by each crop species (cereals, maize, clover and residue from the previous crop field for tillage), respectively, for cereals, maize, clover and tillage fields (% plant cover of each crop species).


**Table A3.** Bird species recorded in the study area (in bold exclusively or partially seed-eaters). \*: The bird species experienced decline in Europe [6,22]; F = farmland specialist, W = primarily woodland species that commonly use farmland [6].

1. For the species listed in Annex II, states that have signed the treaty are required to take the necessary measures for the protection and conservation of these species and their habitats; for the species listed in Annex III, states that have signed the treaty are required to regulate the exploitation of wild fauna and prevent illegal means of capture and killing. 2. I: Species that will be subject of special conservation measures taking into account their habitat to ensure their survival and reproduction in the area of their dispersal; II/1: Species that can be hunted in the geographical sea and land where this Directive applies; II/2: Species that can be hunted only in the Member States, having regard to local laws. 3. 1: Species of global interest for their conservation; 2: Concentrated in Europe and with an unfavorable conservation status; 3: Not concentrated in Europe, but with an unfavorable conservation status; 4: Concentrated in Europe and with a favorable conservation status; w: Category related to winter populations; 4. I: Species with risk of total or at large extent extinction; II: Species that can benefit from the international cooperation for their conservation and management.



**Table A5.** The total number of seeds per m2 estimated for each species of the topsoil seedbank in each studied habitat for winter (left columns) and spring (right columns), respectively.



#### **Table A5.** *Cont.*

Out of the 66 identified species, 35 were commonly detected in both seasons (winter and spring) and are marked with an asterisk (\*); seven species found across all habitats and seasons are additionally shown with the double cross (‡) indicates absence. In **bold** are shown species of (least to highest) significance as food items to rural birds, while significant and highly significant species are besides underlined; classification followed [22,78,79], and field observations.

#### **Appendix B**

List of sources used for seed and plant specimen identification. Seed identification:


http://asis.scri.ac.uk/

• The Ohio State University. Department of Horticulture and Crop Science. Seed IDWorkshop

http://www.oardc.ohio-state.edu/seedid/

• University of Missouri Extension. Missouri Weed Seeds. Department of Agronomy Fred Fishel Kevin Bradley

http://extension.missouri.edu/explore/agguides/pests/ipm1023.htm

• Seeds of Success Collections at the Bend Seed Extractory

http://www.nps.gov/plants/sos/bendcollections/index.htm

• The seed identification web page. Paleoethnobotany Project

http://www.oldthingsforgotten.com/seeds/seeds.htm

• Visual Identification of Small Oilseeds and Weed Seed Contaminants Grain Biology Bulletin No. 3

http://www.grainscanada.gc.ca/Pubs/Grainbio/bulletin3/sows\_03-e.htm Plant specimen identification:


Websites (online databases):


http://extension.usu.edu/weedweb/ident/ID.htm

• University of California, Agriculture and Natural Resources, Statewide IPM Program http://www.ipm.ucdavis.edu/PMG/WEEDS/low\_amaranth.html

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

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

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9575-7