Modelling Impacts of a Municipal Spatial Plan of Land-Use Changes on Surface Water Quality—Example from Goriška Brda in Slovenia

Abstract

Intensive agriculture causes nutrient leaching and accelerates erosion processes, which threatens the good quality status of surface waters, as proposed by the European Union (EU) Water Framework Directive. The purpose of this study was to define the impact of two alternative agricultural land-use change scenarios defined in a Municipal Spatial Plan on surface water quality by using the Agricultural Policy/Environmental eXtender (APEX) model. As experimental area, we chose a small Kožbanjšček stream catchment (1464 ha) situated in the Goriška Brda region in Slovenia. The area, due to favorable conditions for vineyards, is facing increasing deforestation. The change of 66.3 ha of forests to vineyards would increase the sediment, nitrate, and phosphorus loads in the stream by 24.8%, 17.1%, and 10.7%, respectively. With the implementation of vegetative buffer strips as a mitigation measure of the current situation, we could reduce the sediment, nitrate, and phosphorus loads by 17.9%, 11.1%, and 3.1%, respectively, while a combination of the two land-use change scenarios would result in a slight increase of the above-mentioned loads, corresponding to 0.61%, 2.1%, and 6.6%, respectively, compared to the baseline situation. The results confirm that, as we can increase pollution levels with deforestation, we can also reduce water pollution by choosing proper types of land management measures.

1. Introduction

Water bodies in Europe and around the world are subject to many anthropogenic pressures. The food demand increased dramatically in the 20th century due to rapid population growth and industrial development. In response, the traditional extensive agriculture was intensified by mechanized production with increased soil disturbance and inputs of fertilizers and plant protection products. The excessive use of fertilizers led to increased nutrient concentrations in all types of water bodies in specific areas beyond their self-cleaning ability, which led to the degradation of aquatic ecosystems [1,2]. In order to retain their function for humans and natural ecosystems, water bodies, one of the significant components of the biosphere, need to be adequately managed [3]. Due to their ecological and economic importance, the protection of water resources is a high priority for the European Union [4].

The Water Framework Directive (2000/60/EC) adopted for sustainable water management and protection has become a binding and critical document for water management in the European Union (EU) member states [5]. Its main objectives are to achieve a good status of surface and groundwater bodies and to prevent the deterioration of the water status [5].

Agriculture activities can be a source of many emissions that can be displaced by surface and subsurface water flow or leached to water sources, which is often reflected in the deterioration of water quality [6]. In Europe, it is estimated that agriculture contributes between 46% and 87% of total nitrogen and between 20% to 40% of phosphorus inputs in water [7,8]. These figures can significantly vary depending on agriculture type, soil tillage, topography, geology, soil type, and meteorological conditions [9]. Soil erosion and the associated soil particles transport into water resources are among the significant threats from the perspective of sustainable agriculture and surface water quality [10]. When this material is transported along a riverbed by a water stream, it is called sediment load. Floating sediments are also called suspended material. They are mainly clay, silt, and sand particles and are indicators of erosion processes in the catchment area [11]. The occurrence of suspended material in aquatic environments is natural and even has the specific functions of providing substrate and habitat for many organisms as well as of supporting various ecosystem services [12]. On the other hand, excessive soil erosion from agricultural land deteriorates the water status, affecting aquatic organisms by spoiling the water’s chemical and ecological status, causing the loss of useful volume, and changing the hydromorphology and related processes [10,13].

Sustainable agri-environmental measures of agricultural land management can reduce the loss of nutrients and soil erosion, but it is difficult to estimate their costs and impacts on crop production performance and water quality improvements [6]. Appropriate measures should be selected according to their efficiency and the characteristics of the river catchments [14,15]. Best management practices (BMPs) incorporated in the code of good agricultural practices have been implemented worldwide in pursuit of an economically acceptable agricultural activity which would reduce the negative environmental impacts [16]. They cover a range of measures and their combinations (cultivation of soil in parallel with cultivation in greenhouses, permanent greening of agricultural land, crop rotation, conservation tillage, vegetation buffers, terracing slopes, land-use change in a catchment) to protect the aquatic environment while strengthening other ecosystems services, such as soil fertility [14]. These measures mostly include the control of hydrological processes in the catchment, since these processes drive erosion and nutrients transport. Changes in land use/crop type or soil management induce an increase in the soil water-holding capacity as well as in the canopy water storage capacity. These can decrease soil erosion and, consequently, the loss of nutrients into water resources [17].

The evaluation of the impacts of different processes and human activities in the catchments on water quality is well researched [18]. Although field measurements provide information on the actual state of the environment, only a large number of measurements over an extended period can provide relevant conclusions. The evaluation is difficult because of the delayed response of the environment to the implemented measures, the complexity of the hydrological processes, and the variability of the ecosystems [4]. With hydrological models, we can save time and money in the evaluation of the processes, as they can simplify and simulate only processes critical to soil erosion, nutrient leaching, and riverbed sediment load transport. The advantage of a properly calibrated and validated model is its capability to simulate the effects of many combinations of pollutants, soil properties, weather conditions, and agricultural practices in heterogeneous landscapes on water quality and quantity [19,20]. Models have proven to be quite useful in environmental protection planning, also in connection with the implementation of the Water Framework Directive [4]. Models help us to identify critical areas, forecast the effectiveness of measures, and analyze the costs that are relevant to stakeholders, thus providing support in the selection and placement of optimal mitigation measures and agricultural practices in the river catchments [21]. There are several mathematical or catchment models available online that allow modelling the impacts of selected measures at the farm or small river catchment level [3,21]. We chose Agricultural Policy/Environmental eXtender (APEX) [22,23], which is a physical, spatially distributed, and temporally continuous hydrologic model that works on an hourly or daily basis and enables long-term continuous simulations [24].

With this research, we wanted to check different land management scenarios in the area of the Kožbanjšček stream catchment in Goriška Brda, Slovenia. The scenarios tested the effects of deforestation in order to establish new vineyards and the impact of agri-environmental mitigation measures on the quality status of surface waters (soil erosion, sediment load transport, nutrient leaching), as envisaged by the Brda Municipal Spatial Plan. Additionally, the practical utility of the APEX model was tested in humid sub-Mediterranean climatic and flysch soil conditions in Goriška Brda.

2. Materials and Methods

2.1. Case Study Area

The study area covers 1464 ha (14.6 km²) of the Kožbanjšček stream catchment which is part of the larger catchment of the River Reka 3000 ha (30 km²). The catchment is located in the municipality of Brda (72 km²) in the western part of Slovenia on the border with Italy (Figure 1). The topography is diverse and characterized by processes of river erosion. The stream valley is surrounded by numerous hills with an average steepness of the slopes of 33.9% (Figure 1 and Figure 2). The terrain classifies the area as less favored for agriculture production. However, with the cultivation of slopes as terraces, it has been possible to achieve an economically efficient production on slopes with more than 15% inclination. Terraces retain water and reduce soil erosion [25].

The Kožbanjšček stream spring is under the foot of Korada Hill on the northern edge of Brda (Figure 2). It mostly flows through flysch geology and finally discharges into the river Reka in the settlement of Neblo. Neblo is the largest settlement in the area, with a population of about 200 inhabitants, while other settlements are represented by individual smaller hamlets [26].

The soil was formed during the geological period of the Eocene. The flysch consists of alternating layers of sedimentary rocks such as claystone, marl, and sandstone [26]. Flysch is poorly permeable and non-resistant to weathering, which makes the study area vulnerable to erosion and landslides [25]. Eutric brown soil is the dominant soil type, very suitable for agricultural production (Figure 1). Rendzina on limestone is present in the northern part of the catchment. The dominant land uses are forests, vineyards, meadows, orchards, built-up land, and others (

Table 1. Spatial analysis of land use in the Kožbanjščak stream catchment.

ID	Land Use	Area
		km2	ha	%
1100	Arable	0.046	4.64	0.32
1211	Vineyard	1.636	163.60	11.17
1221	Intensive orchards	0.039	3.90	0.27
1222	Extensive orchards	0.346	34.64	2.37
1230	Olive groves	0.039	3.86	0.26
1300	Permanent grassland	1.126	112.64	7.69
1410	Overgrown land	0.234	23.40	1.60
1500	Trees and shrubs	0.175	17.54	1.20
1600	Uncultivated agricultural land	0.033	3.32	0.23
1800	Agro-forestry	0.059	5.92	0.40
2000	Forest	10.522	1052.20	71.85
3000	Urban	0.357	35.73	2.44
6000	Other	0.0005	0.05	0.0003
7000	Water	0.031	3.07	0.21
Total		14.645	1464.50	100

, Figure 1 and Figure 2). Alluvial soils in the valleys formed on carbonate deposits (sand, gravel) are suitable as arable land thanks to their fertility and location; however, vineyards and orchards predominate. The soil surface of all vineyards in the area is covered with grassland to reduce erosion. This is a general cross-compliance measure of the Slovenian Rural Development Programme (RDP) within the EU Common Agricultural Policy (CAP).

The river network of the Kožbanjšček stream and its tributaries is 6996 m in length (Figure 1) with a rain-dominated flow regime. As flysch is impermeable to water, rainfall drains rapidly from the surface, which may lead to the formation of torrential streams. The streamflows are highest in the autumn and winter and lowest in the summer. Groundwater aquifers are almost non-existent due to the impermeable soil [25]. A state-monitoring gauging station is situated at the outflow of the catchment, providing daily information about water discharge.

The area is dominated by warm and sunny sub-Mediterranean climate, which changes on the north-east edge of the catchment (climate divide) into continental climate. In the period 1981–2010, the average annual temperature was 13.1 °C, and the average annual rainfall was 1540 mm (Vedrijan meteorological station). The highest precipitation (140–190 mm/month) is expected in autumn (September–December), defined as wet period, and the lowest (70–100 mm/month) in winter (January–March), defined as dry period [27]. The period between April and July is defined as a transitional period with an undefined weather pattern (100–150 mm/month), although droughts are rare in this period.

The primary sources of water pollution are diffuse and originate from agriculture activity. There are no major point sources in the area, such as municipal or industrial wastewater treatment plants. The most significant problem for a good water quality status (aquatic life) is soil erosion and, consequently, the high concentrations of suspended materials in surface waters, which exceed the limit values (

Table 2. Average and maximum monthly levels of nitrate nitrogen, suspended materials, and total phosphorus (TP) (based on daily observations) at the Neblo outflow station on the Kožbanjšček stream from July 2008 to June 2009.

Month	Nitrate		Suspended Materials		Total Phosphorus
	avg.	max.	avg.	max.	avg.	max.
	(mg/L NO3−)	(mg/L NO3−)	(mg/L)	(mg/L)	(mg/L TP)	(mg/L TP)
January	2.8	6.0	55.5	515.0	0.07	1.06
February	1.9	6.4	56.4	661.0	0.13	1.59
March	1.8	6.0	27.2	135.0	0.62	16.65
April	1.4	12.6	15.3	44.9	0.05	0.24
May	6.6	19.4	31.3	579.0	0.09	1.21
June	2.8	6.2	22.3	451.0	0.06	0.71
July	3.6	10.8	13.3	75.0	0.05	0.27
August	2.6	4.0	16.2	53.9	0.05	0.59
September	2.6	6.5	22.5	147.0	0.04	0.30
October	2.6	7.0	20.8	42.7	0.02	0.10
November	1.5	2.5	43.3	129.0	0.06	0.28
December	2.8	8.3	65.9	162.0	0.08	0.38
Average	2.7	-	32.6	-	0.11	-

Target values in the legislation: suspended materials 25 mg/L; nitrate in drinking water 50 mg/L, nitrate in waters of very good state 14–30 mg/L, nitrate in waters of good state 29–42 mg/L; total phosphorus: 0.2 mg/L in salmonid waters and 0.4 mg/L in cyprinid waters.

).

2.2. The APEX Model

The APEX model is derived from the Environmental Policy Integrated Climate (EPIC) model and is its direct extension. It can model at the farm or at the smaller river catchment level. It was developed to address the environmental problems associated with the intensive livestock production of a project called Livestock and the Environment: A National Pilot Project (NPP) [28]. EPIC and APEX combine the already existing Groundwater Loading Effects of Agricultural Management Systems (GLEAMS), Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS), Carbon and Nutrient Dynamics Ecosystems model (CENTURY), and Agricultural Land Management Alternatives with Numerical Assessment Criteria (ALMANAC) models [27]. The APEX model considers of the following major elements: weather, plant protection products, hydrology, soil erosion/sedimentation, nutrient circulation, plant growth, carbon circulation, soil tillage, soil temperature, and economics (Figure 1) [28,29]. The APEX model enables the simulation of water flow pathways, streamflow, nutrients, sediments, and plant protection products transport between spatially linked sub-areas, such as fields or sub-catchment homogeneous by land characteristics [30,31]. To make its application more user-friendly, interfaces such as SWAPP, ArcAPEX, as well as WinAPEX and iAPEX for Windows and Mac were developed [22]. In this study, we used ArcAPEX 10.3.2 version released in May 2017 that runs in ArcGIS 10.3.1 (Geographic Information System) and is compatible with APEX v.1501.

2.3. Database

The APEX model requires a wide range of input data to run. These include weather data, soil properties by horizons, land management and general catchment characteristics (

Table 3. Database set-up.

Data			Vir
Meteorological Stations	Brdice, Vedrijan	Daily precipitation	ARSO
	Bilje	Daily temperature (min., max.)
		Sun hours
		Relative air humidity
		Wind
Spatial data		Digital elevation model	GURS
		Land use	MKGP
		Soil type	MKGP
		Soil horizons properties	BF
		River network	APEX
		Gauging stations	ARSO
Land use management		Crop rotations	BF
		Fertilization	BF
		Land management	BF
		Forest	ZGS
Calibration (Neblo station)		Daily streamflow discharge	ARSO
Validation (Neblo station)		Water quality (nitrate, phosphorus, suspended materials)	BF

Note: ARSO—Slovenian Environmental Agency, GURS—Surveying and Mapping Authority of the Republic of Slovenia, MAFF—Ministry of Agriculture, Forestry, and Food, BF—Biotechnical Faculty, ZGS—Slovenian Forest Service.

). All daily weather data for the simulation period were obtained from the Slovenian Environmental Agency web-based database [27]. For daily rainfall, we used data from the nearest precipitation stations of Brdice and Vedrijan. Missing daily weather data at the Brdice station were replaced by the Vedrijan station data. For all other weather variables, we used data from the Bilje meteorological station.

2.4. Model Set-Up

To divide the catchment into sub-areas and present the river network distribution, we used a 25 × 25 m digital elevation model (DMV) spatial layer. The threshold for the minimum area required to identify sub-areas and determine their number and size [16] was set at 58 ha and 928 cells to obtain adequate accuracy of the river network. The 1462 ha-large study area was divided into six sub-areas that are hydrologically interconnected in terms of surface runoff, pollutant transport, subsurface flow, and other processes (Figure 1). Since the model works by assigning only one land-use class, soil type, and slope to each subarea, we decided to define the subareas on the basis of actual land use and soil maps. Land use in sub-areas 1–5 was set as forest, and that in sub-area 6 as vineyards. Daily agricultural task calendars for the selected land uses in the area were prepared and included in the model (

Table 4. Agricultural management task calendar for vineyard (GRAP) and forest (FRSD) land uses in the APEX model.

Crop	Month	Day	Operation	Operation Type	Fertilizer Application (kg N, P/ha)	HE	ORHI
GRAP	1	10	Harvest without kill	Cut of shoots		0.01	0.33
	3	4	Plant in rows	4000 vines/ha
	4	1	Fertilizer	Elemental N	20
	4	1	Fertilizer	Elemental P	23
	5	1	Fertilizer	ORG-N (mulch)	144	1	0.9
	5	1	Fertilizer	ORG-P (mulch)	11
	6	14	Fertilizer	Elemental N	50
	6	15	Fertilizer	ORG-N (mulch)	144	1	0.9
	6	15	Fertilizer	ORG-P (mulch)	11
	7	15	Fertilizer	ORG-N (mulch)	144	1	0.9
	7	15	Fertilizer	ORG-P (mulch)	11
	8	20	Fertilizer	ORG-N (mulch)	144	1	0.9
	8	20	Fertilizer	ORG-P (mulch)	11
	9	5	Harvest without kill	Harvest fruit		1	0.2
	11	20	Fertilizer	ORG-N (mulch)	346
	11	20	Fertilizer	ORG-P (mulch)	41
FRSD	1	1	Plant in rows	Transplanter, 100 t/ha
	10	31	Plough, cultivate, other	Thin

Key: FRSD—broadleaf forest, GRAP—vineyard, HE—harvest efficiency, ORHI—override harvest index, ORG-N—organic nitrogen, ORG-P—organic phosphorus.

). The outflow points were manually determined to obtain the most homogeneous sub-areas. The official state hydrological and water quality monitoring station Neblo (daily flow discharge monitoring) was determined as the main outflow point (streamflow discharge) and calibration point (streamflow, sediment load, and nutrients loads).

To present results consistent with the APEX terminology, we show results of river water discharge as daily or monthly streamflow at the catchment outflow, suspended materials as sediment load transported out of the catchment, nitrate as nitrate nitrogen load transported out of the catchment, and phosphorus as total phosphors load transported out of the catchment.

The simulation period from 1 April 1991 to 1 September 2009 was divided into warm-up, calibration, and validation periods, with daily time-step modelling. Since water is an essential driving force for river catchment processes, proper modelling of hydrological conditions (water balance) was crucial for nutrients and sediment load calibration [32,33]. The model was calibrated and validated for daily streamflow discharge, sediment load, nitrate and total phosphorus loads. For the hydrological calibration, warm-up (1991–1997), calibration (1998–2005), and validation (2006–2008) periods were defined. The calibration of sediment and nutrients loads was performed on daily measured data obtained during a one-year monitoring campaign (1 July 2008 to 30 June 2009) as in one of the previous studies [34]. Validation was not carried out, due to lack of data. The performance of the model was evaluated on daily and monthly time scales.

Sensitivity analysis, calibration, and validation of the model were performed automatically with the open-source Agricultural Policy Environmental eXtender–Auto-Calibration and UncerTainty Estimator (APEX-CUTE) tool. The tool was developed primarily for the calibration of the APEX model and uses the Dynamically dimensioned search (DDS) algorithm for parameter optimization [35,36]. The sensitivity analysis was done for the sub-area 6 outflow, setting the number of intervals to 100, i.e., 100 iterations per parameter and another 100 per parameter combinations. The tool allows the analysis of 160 parameters describing different processes in the catchment. In the analysis, the parameters for sensitivity analysis were selected based on the literature [29,36]. Those identified by the program as sensitive and used in the calibration process are presented in the results (

Table 5. Parameters identified in the sensitivity analysis used for the calibration and validation of streamflow, sediment load, nitrate nitrogen and total phosphorus loads.

Parameter		Description	Default Value	Calibrated Value
Streamflow	APM	Peak runoff rate—rainfall energy adjustment factor	1	0.559
	FCW	Floodplain width/channel width in m/m	1	30.465
	QCF	Exponent in watershed area flow rate equation	0.5	0.425
	PARM12	Soil evaporation coefficient	1.5	2.145
	PARM17	Soil evaporation–plant cover factor	0.1	0.376
	PARM20	Runoff curve number initial abstraction	0.2	0.156
	PARM40	Groundwater storage threshold	0.25	0.198
	PARM42	SCS curve number index coefficient	0.5	0.809
	PARM46	RUSLE C-factor coefficient	0.5	0.754
	PARM49	Maximum rainfall interception by plant canopy	7.5	14.402
Sediment load	APM	Peak runoff rate—rainfall energy adjustment factor	1	0.275
	RCC0	USLE Crop Management Channel Factor	0	0.057
	RCHK	USLE erodibility channel (K) factor	0.300	0.0001
	UPN	Manning’s N for upland	0.150	0.100
	PARM18	Sediment routing exponent of water velocity function for estimating potential sediment concentration	1.5	1.00
	PARM19	Potential sediment concentration when flow velocity is 1.0 m/s	0.05	0.01
	PARM46	RUSLE C factor coefficient in exponential residue function in residue factor	0.5	0.962
	PARM47	RUSLE C factor coefficient in exponential crop height function in biomass factor	1	0.581
Nutrients load	PARM14	Nitrate leaching ratio	0.2	0.849
	PARM29	Biological mixing efficiency	0.1	0.185
	PARM70	Microbial decay rate coefficient	1	1.069
	PARM72	Volatilization/nitrification partitioning coefficient	0.15	0.433
	PARM74	Partitions Nitrogen flow from groundwater	10	17.179

SCS—Soil Conservation Service, RUSLE—Revised Universal Soil Loss Equation.

). The calibration process was automatic, setting the maximum number of iterations to 1000. Text files with daily time series of measured data were prepared, according to which APEX-CUTE calibrated the simulated values. APEX-CUTE determined the optimal parameter values according to the best combination of Nash–Sutcliff simulation efficiency (E_NS) and percent bias (PBIAS) after altering the parameter values and after 1000 repetitions of autocalibration [37,38]. The best set of final parameter values was then used to test the base model objective function in the validation period.

2.5. Objective Functions

Several statistical and graphical methods are available to evaluate the predictive ability of the model. It is important to use several indicators at one time in the assessment, as they reveal different uncertainties of the model [36]. Generally recognized objective functions were used for model calibration and validation performance tests, such as E_NS and PBIAS, used also by APEX-CUTE, as well as others (coefficient of determination–R², effective root-mean-square error–RMSE, and standard deviation) [19].

R² (unitless) describes the portion of the total variance in the measured data that can be explained by the model. The range is from 0.0 (poor model) to 1.0 (perfect model). A value of 0 means that the model replicates none of the variances in the measured data, and a value of 1 means that the model predictions replicate all the variances in the measured data. The value acceptable for hydrological studies is 0.5 (streamflow). E_NS (unitless) is used to evaluate the performance of the hydrological model. It measures how well the simulated results predict the measured data. Values range from negative infinity (poor model) to 1.0 (perfect model). A value of 0.0 means that the model predictions are as accurate as those made by using the measured data average. A value higher than 0.0 means that the model is a better predictor of the measured data than the measured data average. A significant disadvantage is that the differences between the measured and the simulated values are calculated as squared values, and this places emphasis on peak flows. As a result, the impact of larger values in a time series is overestimated, whereas lower values are neglected. Values should be above zero to indicate a minimally acceptable performance. The value acceptable for hydrological studies is 0.5 (streamflow). RMSE is determined by calculating the standard deviation of the points from their true positions, summing up the measurements and then taking the square root of the sum. RMSE is used to measure the difference between the flow values simulated by a model and the actual measured flow values. Smaller values indicate a better model performance. The range is between 0 (optimal) and infinity. The value acceptable for hydrological studies is 0.7 (streamflow). PBIAS (%) measures the average tendency of the simulated flows to be larger or smaller than their observed counterparts [19]. The optimal value is 0, positive values indicate a model bias toward underestimation, and negative values indicate the opposite. The value acceptable for hydrological studies is 25% (streamflow).

2.6. Scenarios

We developed two alternative land management scenarios to assess their impact on surface water quality. The scenarios were introduced in the model by changing the parameters in the plant database, tillage, calendar of agriculture management tasks, and properties of the sub-areas. The simulation period for the alternative scenarios was 1 April 1991 to 1 September 2009, excluding the warm-up period (1991–1994) from the analysis. The impact of different agricultural land management methods was evaluated by comparing the simulation results of the baseline scenario with those of the alternative scenarios on an average annual basis. All scenarios included terraces as a base land management practice in the area. The comparison of the results between the scenarios is presented as a percentage change in the average annual sediment, nitrate nitrogen load, or total phosphorus load in the streamflow (%).

2.6.1. Land Use Change (LUC5)

Part of the Municipal Spatial Plan of the Brda municipality comprises initiatives for deforestation and conversion to vineyards, which would preserve the characteristic cultural landscape and promote tourism [39]. This process could accelerate soil erosion and landslides and affect the water quality status of surface watercourses and aquatic life due to the removal of natural vegetation [39,40].

In the LUC5 scenario, we changed forest (FRSD) land use in the sub-area 5 to vineyards (GRAP), on the basis of topographic (less than 50% slopes) and soil properties (eutric brown soil on flysch) most suitable for the establishment of new agricultural land. Land-use change increased the total area of the vineyards in the model from 188.2 ha to 254.5 ha, representing 17% of the area, while the forested area in the model was 83% of the Kožbanjšček stream catchment. The parameters and technology of vineyard management were the same as in the baseline scenario, including terracing.

2.6.2. Vegetative Buffer Strips (VBS6, VBS56)

Since flysch soils are prone to erosion, measures should be taken to reduce its effects in areas of newly established agricultural land. Vegetative buffer strips are a simple and cost-effective measure to reduce the transfer of eroded soil particles and nutrients to water sources. With this measure, a strip of natural or sown vegetation is established between the source of pollution and the watercourse, catching pollutants present in the surface runoff and increasing the infiltration of water into the soil [41].

The latest version of APEX offers us two options for simulating vegetative buffer strips. The first approach requires a precise spatial definition of the buffer zones and their hydrological relationship with the surfaces that are the source of the eroded material. The model requires a definition of strip width, vegetation type, soil conditions, topographic features (slope, roughness), proportion of concentrated stream flowing through the strip, vegetation growth, and buffer strip management (cutting, fertilization, irrigation) [42]. In the second, more straightforward approach, chosen for this scenario, we created virtual buffer strips, which can be placed in any sub-area without knowing their exact locations. This method is suitable for evaluating the effectiveness of the strips at the catchment level. The physical characteristics of the strips such as slope, soil, and vegetation type cannot be altered, but we can determine the concerned sub-areas, as well as the proportion of sub-areas from which the surface runoff flows through the strip (parameter BCOF) and its width (parameter BFFL) [43]. We placed vegetative buffer strips on all vineyards (GRAP) in sub-area 6, which is also the most critical in terms of soil erosion and nutrient runoff (VBS6 sub-scenario). We assumed that 20% of all surface runoff from the subarea flowed through 5 m-wide vegetative buffer strips [44]. Additionally, sub-scenario VBS56 was prepared to evaluate the impact of vegetative buffer strips on soil erosion reduction in both sub-areas 5 (vineyards) and 6 (66 ha of the forest changed to vineyards). The parameters and technology of vineyard management were the same as in the baseline scenario, including terracing.

3. Results and Discussion

3.1. Sensitivity Analysis, Calibration, and Validation

3.1.1. Streamflow

In the process of streamflow sensitivity analysis, we identified 10 sensitive parameters (Table 5). After autocalibration, the agreement between simulated and measured values (

Table 6. Comparison of mean, standard deviation, minimum and maximum daily mean streamflow (m³/s) between measured and simulated values for the baseline model in the calibration (1998–2005) and validation (2006–2008) periods at the Kožbanjšček catchment main outflow.

Statistics	Calibration Period		Validation Period
	Simulation	Measurement	Simulation	Measurement
Average	0.32	0.32	0.16	0.14
Standard deviation	1.00	0.82	0.42	0.36
Minimum	0.00	0.00	0.00	0.00
Maximum	26.63	13.40	7.55	4.88

and

Table 7. Statistical evaluation of simulated mean daily and monthly streamflows for the calibration (1998–2005) and validation (2006–2008) periods at the Kožbanjšček catchment main outflow.

Objective Function	Daily Time-Step		Monthly Time-Step
	Calibration	Validation	Calibration	Validation
R2	0.498	0.268	0.809	0.806
ENS	0.233	−0.145	0.705	0.761
RMSE	4.226	2.272	0.186	0.074
PBIAS (%)	0.351	−13.708	0.805	−13.878

) of mean daily streamflow proved to be acceptable [19]. The PBIAS indicated less than 1% deviation from the measured values and was within the recommended range for both daily and monthly calibration. Slightly positive values indicate an underestimation of the flow. The coefficient of determination (R²) values also was satisfactory [19]. The E_NS values for the daily calibration were slightly lower and did not reach the satisfactory value of 0.5 [19]. The reason may lie in the sensitivity of the E_NS to extreme values [19]. Calibration on a monthly base returned values of 0.81 (R²), 0.71 (E_NS), 0.19 (RMSE), and 0.81 (PBIAS), indicating good and very good agreement between the measured and the simulated data.

The validation results at the daily and monthly levels returned lower values, as found by other authors [29]. For daily validation, R² and E_NS values were outside and PBIAS was within the recommended range (Table 7). For monthly validation, the values of 0.81 (R²) and 0.76 (E_NS) showed a good match between simulated and measured values. PBIAS was within the range of good values. The peaks coincided, with occasional underestimation of the flow in the calibration period and overestimation of the flow in the validation period (Table 7).

Satisfactory calibration and validation results are essential for further work with the model. Hydrological calibration and validation were not ideal. One of the reasons may be that we used in the calibration only a limited number of the parameters primarily recommended by the sensitivity analysis, while the latest version of APEX-CUTE offers 160 different parameters for calibration (Table 5). The parameters were also optimized at the whole catchment level, while in reality, there were differences between individual sub-areas that would require individual calibration. Figure 3 indicates the impact of the model parameters, set up in wetter hydrological conditions of the calibration period, on the modelling results of the validation period when hydrological conditions were clearly drier (lower average monthly streamflow). Parameters calibrated in wetter or dryer periods showed a lower agreement with those determined in dryer or wetter periods, respectively. Last but not least, the plants had an impact on the water balance in the catchment; however, the model was not calibrated for biomass and crop production as recommended by the literature [29] due to unavailability of data.

3.1.2. Sediment Loads, Nitrate Nitrogen, and Total Phosphorus Loads

In the process of the sensitivity analysis, we identified eight sensitive parameters for the sediment load and five sensitive parameters for the nutrients load (

Table 8. Comparison of the average daily value, standard deviation, minimum and maximum sediment loads, total phosphorus, and nitrate nitrogen between measured and simulated values after calibration of the base model (2008–2009) at the Kožbanjšček stream main outflow.

Statistics	Simulation	Measurement
	Sediment Load (t/day)
Average	3.01	2.93
Standard deviation	19.19	26.33
Minimum	0.00	0.00
Maximum	263.5	407.32
	Total phosphorus load (kg/day)
Average	6.30	6.40
Standard deviation	36.02	47.03
Minimum	0.00	0.00
Maximum	449.39	638.52
	Nitrate-Nitrogen load (kg NO3-N/day)
Average	18.66	34.08
Standard deviation	82.17	102.86
Minimum	0.00	0.01
Maximum	922.60	851.12

and

Table 9. Statistical evaluation of simulated daily sediment load and total phosphorus and nitrate nitrogen loads for the calibration period (2008–2009) at the Kožbanjšček stream main outflow.

Objective Function	Calibration Period
	Sediment	Total Phosphorus	Nitrate Nitrogen
R2	0.277	0.158	0.543
ENS	0.236	0.023	0.520
RMSE	22.851	46.491	71.492
PBIAS (%)	−2.770	1.516	45.245

). In the calibration of the model, we used measured daily values for the period between 1 July 2008 and 30 June 2009. As autocalibration with APEX-CUTE did not return satisfactory results, the parameters were changed manually to better match the simulated and measured data.

The PBIAS results for sediment load indicated that the predictive ability of the model was very good (Table 8 and Table 9) [19]. Even the average value did not deviate significantly from the measured data (Table 8). The R² and E_NS results indicated lower agreement, but we have to keep in mind that the calibration period was short, and the E_NS values improved with temporally denser data. A visual comparison showed that the simulated values followed seasonal fluctuations, but over time the peak values underestimated or overestimated the sediment load in the water (Figure 4). During the one-year water sampling campaign, earthworks for the construction of new cultivation terraces were performed in the catchment and could not be captured in the model set-up. As mentioned in the scenarios’ description, terraces were included in modelling as a base practice in the area.

Total phosphorus in water represents both dissolved and particle-bound phosphorus. The APEX provides separate results for phosphorus as a mineral and as organic phosphorus, which were for this study summed, to obtain a value for total phosphorus. The calibration results of total phosphorus indicated a very good agreement between the simulated and the measured values in the case of calculated PBIAS (1.5), while the R², E_NS, and RMSE values were lower (Table 8 and Table 9, Figure 4). These results are closely related to the fact that phosphorus is mostly transported bound to soil particles in surface run-off as part of the soil erosion process [45,46].

The APEX-CUTE tool does not offer the possibility of directly calibrating nitrate nitrogen, and only total mineral nitrogen can be determined. In normal conditions in the area, nitrite and ammonium nitrogen in the river water account for less than 10% of the total mineral nitrogen. For the calibration, we had to assume the approximation that the nitrate nitrogen content, for which we had one-year daily measured data, was equal to the total mineral nitrogen content. The APEX model expresses nitrogen content in water as nitrate nitrogen (NO₃-N). This approximation in the calibration phase led to a slight underestimation of less than 10%, calling for an upgrade of APEX-CUTE, allowing calibration with nitrate nitrogen.

The objective functions of E_NS (0.52) and R² (0.54) were in the satisfactory range. PBIAS (45.25) was still in the acceptable range (Table 9) [19], although the average value and standard deviation of the simulated data differed from the measured values (Table 8). The reason for the larger deviation observed could be the low average concentrations of nitrate in water (2.7 mg/L), for which a small difference of 0.5–1 mg/L represented a 15%–30% deviation, which is not harmful to the environment in this area. A visual comparison showed that the model indicated seasonal changes of nitrate in water but, in general, underestimated the values.

When the model was calibrated for the sediment load and nutrients load, we faced some challenges in achieving satisfactory results. In addition to earthworks, one of the reasons could be the way APEX attributed land use and other properties to sub-areas. In the test area, APEX attributed only two dominant land uses to the subareas, forest (sub-areas 1–5) and vineyard (sub-area 6), which was far from the actual heterogeneous land use situation. The model excluded from the simulation the orchards, meadows, olive groves, fields, and urban areas; thus, it was impossible to completely simulate the actual flow of nutrients and soil particles (erosion) into the stream. Besides, we only had a one-year set of daily measured data available for calibration. The reasons for the poor calibration results may also lie in the measurement uncertainty of the measured data, the deficient time series of the input data, and the limited modelers’ experience in parameter tuning [47].

3.2. Scenarios Results

The Land Use Change (LUC5) scenario would have a significant impact on the water quality status in the Kožbanjšček stream (

Table 10. Impact of alternative scenarios on the change (%) in average annual sediment load (t/year), total phosphorus load (kg/year), and nitrate nitrogen load (kg/year) at the Kožbanjšček stream catchment main outflow in comparison to the baseline scenario.

Parameter (Variable)	Base Scenario Load	Scenario–Percentage Change (%)
		LUC5 Sub-Area 5	VBS6 Sub-Area 6	VBS56 Sub-Area 5 and 6
Sediment (t/year)	342.2	+24.8	−17.9	+0.62
Nitrate nitrogen (kg/year)	21464.9	+17.1	−11.1	+2.1
Total phosphorus (kg/year)	1679.8	+10.7	−3.1	+6.59

Note: LUC5—scenario of modified land use in sub-area 5 from forest to vineyard; VBS6—scenario of vegetation buffer zones in sub-area 6; VBS56—scenario of modified land use in sub-area 6 from forest to vineyard and placement of vegetative buffer zones in sub-areas 5 and 6.

). The land-use change (66.3 ha) in sub-area 5 from forest to vineyard would increase the average load of sediment and nitrate nitrogen and total phosphorus loads in the streamflow by +24.8%, +11.1%, and +10.7% respectively, in comparison to the baseline scenario. Similar impacts of agricultural intensification on slopes in Mediterranean catchments were confirmed by many authors [48,49,50]. Increased nitrate and phosphorus levels in stream water reflect changes in agricultural management consisting in the fertilization of vineyard soils, while increased sediment loads reflect increased soil erosion due to the construction and maintenance of cultivation terraces and the removal of natural vegetation (forest), which has a strong protective function against erosion processes.

On the other hand, the implementation of vegetative buffer strips (VBS) would decrease the average annual loads of sediment, nitrate nitrogen, and total phosphorus in the streamflow by −17.9%, −11.1%, and −3.1%, respectively. The differences between the average annual values for all variables were statistically significant (95% confidence interval).

The scenario VBS56 showed good results and represents a viable compromise, since the increase of sediment load (+0.62), total phosphorus load (+6.59%), and nitrate nitrogen load (+2.1%) was minimal compared to that of the LUC5 scenarios.

The scenarios results incorporate some uncertainties which need to be addressed. Since the exact spatial location of the buffer strips in the model was unknown, it is challenging to define the actual performance of the strips when placed in nature. The APEX model simplifies by simulating a steady flow across the strip and does not account for surface differences [44]. In reality, its performance is strongly influenced by the slope, surface terrain type, and soil properties (erodibility) at micro-locations. Adjustments are essential for rugged terrain, and the holding capacity of the buffer strips is not evenly distributed due to the surface runoff often concentrated at one point. Considering this, it is possible to adjust the width of the strips according to the given area characteristics and ensure their optimal performance. In the case of steep slopes, the width of the strips should be increased accordingly [44]. The most effective buffer strips are a mixture of tree, shrub, and herbaceous plants, while the APEX model simulates strips only using Bermuda grass properties.

Nutrients and sediment loads in waters depend on a wide range of factors, such as weather factors, hydrological processes in the catchment, land use, soil properties, and slope [50]. The time and quantities of fertilizer application, the time and practice of soil cultivation, and the time of crop planting have a significant influence on the growth of the vegetative cover protecting the soil surface [51,52]. Due to the nature of the model operation, particular adjustments were made for the catchment in terms of heterogeneity of land use, soil, and slope, which could significantly affect the model’s results and calibration [53,54]. In reality, cultivation techniques vary from parcel to parcel and from year to year, while in the simulation, they were constant throughout the years of the simulation period and for all lands with the same vegetation cover. The subdivision of the catchment into even smaller sub-areas might better represent the variability in land use, but this is not always possible.

It is worth mentioning that, in comparison to the majority of the EU member states, the Slovenian Ministry for Agriculture, Forestry, and Food, as the head policymaker, allows the reduction of forest land in favor of agricultural land. The main reason for this policy is that forest covers more than 58% of Slovenia (EU 28–8%), and agricultural land represents 23% of the land surface (EU 28–43%), of which cropland represents only 9% (EU 28–22%). Although the Slovenian Rural Development Programme (RDP) financially supports the sustainable development of forests, the diverse topography with steep slopes, high altitude, and unfavorable geology (karst) is the best protector of forests in Slovenia, Landscape is also protected by wetland flood plains, protected forests and Natura 2000 (37% of Slovenia surface). In addition, less-favored areas, as defined by CAP, cover more than 72% of Slovenia.

3.3. Integrated Spatial Planning and Environmental Assessment

The most important economic activities in the Goriška Brda area are closely linked to agriculture and tourism. Upper Brda (including Kožbanjšček catchment) is under demographic threat because of challenging conditions for farming (steeper slopes), poor transport connections, and lack of jobs. This combination has resulted in the abandonment of agriculture, emigration of the population. and overgrowing of agricultural land. Therefore, the development interests of the municipality include the promotion of viticulture, fruit growing, and, in relation to them, tourism. The Municipal Spatial Plan envisages deforestation and the development of new vineyards in order to preserve the agricultural activity, maintain the population in villages, and protect the characteristics of the cultural landscape [38]. As indicated by our results, the implementation of this plan would have a significant impact on the environment.

The scenarios we proposed were evaluated according to the European and Slovenian legislation, which defines the limits and recommends acceptable levels of chemicals and materials in surface waters (

Table 11. Comparison of the effects of land-use scenarios on the average annual concentration of sediment, nitrate nitrogen, and total phosphorus (mg L⁻¹) at the outflow of the river Kožbanjšček—values were recalculated from the average daily loads using daily streamflow.

	Average Annual Observed Concentration (mg L−1)
	Sediment (guideline) b	Nitrate-Nitrogen (NO3-N) a,c			Total phosphorus b
		Drinking Water	Very Good State	Good State	Salmonid Waters	Cyprinid Waters
Limit and guide concentrations	25	11.36	3.2–7.0	6.5–9.5	0.2	0.4
Measured concentrations	32.6	0.61			0.109
Scenarios d
LUC6	40.7	0.71			0.120
VBS5	26.8	0.54			0.106
VBS56	32.8	0.62			0.116

^a 1 mg L⁻¹ of NO₃-N (nitrate nitrogen) is equivalent to 4.4 mg L⁻¹ of NO₃⁻ (nitrate); ^b sediment and total phosphorus (official gazette of Slovenia No. 46/2002; ^c directive 2006/44/EC), nitrate nitrogen (Water Framework Directive 2000/60/EC, official gazette of Slovenia No. 14/2009); ^d LUC5—scenario of modified land use in sub-area 5 from forest to vineyard, VBS6—scenario of vegetation buffer zones in sub-area 6, VBS56—scenario of modified land use in sub-area 6 from forest to vineyard and placement of vegetative buffer zones in sub-areas 5 and 6.

). The measured content of suspended materials (25 mg/L) exceeded, annually and at times monthly, the value established in the decree on the quality required for surface waters to support freshwater fish life [55] (Table 2). The measured total phosphorus content also exceeded the limit values for salmonid (0.2 mg/L) and cyprinid (0.4 mg/L) water in individual cases, whereas the measured nitrate content never exceeded the limit value indicated in the Water Framework Directive (WFD) [5] (50 mg NO₃⁻/L). Table 11 shows the modelled average annual content of substances in water after the implementation of the LUC5 and VBS scenarios. The results indicate that the implementation of the LUC5 scenario would significantly increase erosion processes and, thus, the amount of sediment in surface water. It would also affect nutrients, which, however, would remain within the accepted limits on an annual basis. The results showed that the implementation of vegetative buffer strips is necessary in the areas of vineyards expansion in order to mitigate the effects of the more intensive use of soil. Although the Kožbanjšček stream is not classified in the aforementioned decree as a habitat important for the life of freshwater fish species, it is home for other important and endangered animal species, such as the dragonfly Cordulegaster heros and the river crab Austropotamobius pallipes, which could be adversely affected by the deterioration of the water status, in disagreement with the current water and nature protection policy.

If the plan is to be implemented, new vineyard areas should be appropriately spatially positioned and implemented with selected mitigation measures in order to reduce the impact of soil erosion on water quality. The construction of terraces is planned for all slopes [38], as it can reduce erosion by more than a third [56,57]. The downside of terraces is their expensive construction and maintenance, but in this area, they have been part of the cultural landscape for centuries. The plan also requires the installation of buffer strips in the vicinity of watercourses. This measure is less expensive and less maintenance-intensive, with favorable side effects such as increased biodiversity of natural habitats and reduction of water temperature due to shading, which increases the dissolved oxygen content required by many aquatic organisms. If the strips are large enough, they can provide corridors for wildlife movement [44]. On the other hand, the installation of buffer strips means loss of productive agricultural land. Improved utilization of space could be achieved by implementing the most effective agri-environmental measures or supporting good agricultural practices in the most problematic locations in the catchment.

3.4. APEX Model Assessment

The ArcAPEX is a user-friendly graphical interface that makes it easy for the user to work with the model. The suitability of the APEX model depends mostly on the user’s prior knowledge and experience. A large number of parameters and inputs that can be modified in the process of model building can cause problems to an inexperienced user [20] but allows a more accurate modelling of processes in a catchment. One of the advantages of the APEX model is its flexibility in designing production technologies. It can simulate mulching, mixed crops, conservation tillage, and many other agricultural practices that are not possible to consider with comparable tools. The APEX-CUTE tool shortens the calibration time of the model, as it enables simultaneous calibration and validation of the model and calculates model performance statistics.

We noticed some drawbacks when using the APEX model and the APEX-CUTE tool. Sub-areas in APEX are defined by uniform land use, soil type, slope, and cultivation technologies, which reduces modelling accuracy. This was the biggest drawback of the model when used, in all situations, in the heterogeneous Kožbanjšček stream catchment. The model is thus more suitable for modelling farms and smaller river catchments, which are homogeneous in most of their characteristics. For the simulation of heterogeneous catchments, it is more appropriate to integrate this model with the SWAT model [28,30], which subdivides the subareas (sub-catchments) into individual hydrologically response units (HRU). Another limit was the restricted plant database that does not contain all the agricultural plants (olive, peach, and cherry) that are common in the Mediterranean area as well as in the modelled area. Simulation with trees or permanent crops in the adult stage and in conditions of full of fertility (forest, vineyard) was not possible at one click. The problem can be partially solved by setting the Seeding Rate (SDW) parameter, which determines the initial plant biomass at the start of the simulation and the parameter RTN0, which indicates the number of years of soil tillage before the simulation is started. The latter affects the amount of available carbon and nitrogen in the soil. The most significant drawbacks faced by modelers are the lack of proper instructions for using the latest versions of the model, as well as several bugs that complicate the user experience [37].

4. Conclusions

With this research, we wanted to test different land management scenarios in the area of the Kožbanjšček catchment area in Goriška Brda, Slovenia. The scenarios examined the effects of deforestation while establishing new vineyards on the environmental pollution of surface waters (erosion, sediment load, nutrients loads). At the same time, we tested the APEX computer model in selected Mediterranean climatic and flysch soil conditions.

We found that the change of the forest land use into vineyards would, in the case of the Kožbanjšček stream, increase sediment load and nutrients loads in the water. At the same time, the placement of vegetative buffer strips would almost minimize the impact of soil erosion on surface water quality. Spatial planning at the local level, where private or public interests meet, has, therefore, to involve all stakeholders in the area (landowners, environmental non-governmental organization (NGOs), local and national policymakers) to minimize conflicting situations affecting the environment. We confirmed that APEX is suitable for use in the climatic and soil conditions of Slovenia. However, small and, in terms of land characteristics, homogeneous sub-areas need to be identified when applying the model.

The results obtained contain some degree of uncertainty. They are affected by the computational capability of the model, the modeler’s experience in integrated catchment modelling, as well as the uncertainty in the measured data and calibration results. A good representation of the research area in terms of land use, soil type, and agricultural management is key to an accurate model prediction. One of the uncertainties was due to APEX-CUTE allowing calibration with respect to total mineral nitrogen, while only nitrate nitrogen content data were available.

Future work should involve the integration with the SWAT model and a broader range of different agri-environmental measures contributing to the reduction of factors negatively affecting the quality of surface water resources. An economic analysis of the implementation of cost-effective agri-environmental mitigation measures in agricultural landscapes would even be more critical and crucial for decision-makers.