How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes

Barton, C. Michael; Ullah, Isaac; Heimsath, Arjun

doi:10.3390/land4030578

Open AccessArticle

How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes

by

C. Michael Barton

^1,*,

Isaac Ullah

^1,† and

Arjun Heimsath

^2,†

¹

School of Human Evolution & Social Change and Center for Social Dynamics & Complexity, PO Box 872402 SHESC, Arizona State University, Tempe, AZ 85287, USA

²

School of Earth and Space Exploration (SESE), Arizona State University, PO Box 876004, Tempe, AZ 85287, USA

^*

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Land 2015, 4(3), 578-606; https://doi.org/10.3390/land4030578

Submission received: 23 April 2015 / Revised: 24 June 2015 / Accepted: 1 July 2015 / Published: 14 July 2015

(This article belongs to the Special Issue Agent-Based Modelling and Landscape Change)

Download

Browse Figures

Versions Notes

Abstract

:

We use the hybrid modeling laboratory of the Mediterranean Landscape Dynamics (MedLanD) Project to simulate barranco incision in eastern Spain under different scenarios of natural and human environmental change. We carry out a series of modeling experiments set in the Rio Penaguila valley of northern Alicante Province. The MedLanD Modeling Laboratory (MML) is able to realistically simulate gullying and incision in a multi-dimensional, spatially explicit virtual landscape. We first compare erosion modeled in wooded and denuded landscapes in the absence of human land-use. We then introduce simulated small-holder (e.g., prehistoric Neolithic) farmer/herders in six experiments, by varying community size (small, medium, large) and land management strategy (satisficing and maximizing). We compare the amount and location of erosion under natural and anthropogenic conditions. Natural (e.g., climatically induced) land-cover change produces a distinctly different signature of landscape evolution than does land-cover change produced by agropastoral land-use. Human land-use induces increased coupling between hillslopes and channels, resulting in increased downstream incision.

Keywords:

landscape; land-use; modeling; Spain; erosion; surface process; archaeology; social-ecological system; Mediterranean

1. Introduction

Characteristic features of many Mediterranean landscapes are deeply incised, intermittent watercourses, termed barrancos in Spanish. These can range from modest gullies a meter deep to extensive drainage systems that extend over tens of kilometers and tens of meters deep. While some deeper barrancos may intersect local water tables, most only carry water periodically during or after significant precipitation. Some barrancos are old and are related to characteristics of underlying lithologies, bedrock structure, long-term climate-driven changes in vegetation cover, and regional drainage networks [1]. However, many are clearly much younger, and are incised into unconsolidated sediments or soft calcareous bedrock [2,3,4]. Many older bedrock barrancos also show evidence of recent incision [1,5,6].

There is a widespread consensus that anthropogenic factors—especially agropastoral land-use—played a significant role in Holocene erosion and soil loss throughout the Mediterranean, although there remains considerable debate over the relative causal importance of human and natural processes at different temporal and spatial scales [6,7,8,9,10]. Certainly, the Mediterranean region today and in the recent past is characterized by high sediment transport levels, as a result of both sheet erosion and incision [4,6]. There is also evidence of significant episodes of erosion, including incision, at various times in the historic and prehistoric past that seem coeval with changes in agropastoral land-use patterns [7,11,12,13].

While there has been considerable study of the impacts of land-use on and hillslope and rill erosion in the Mediterranean, the relationships between land-use and gullying and barranco formation are less well understood [4,6,14]. Moreover, quantitative studies of the effects of farming and herding practices on gully incision have been largely empirical and limited to short-term processes (from single events to several decades) that are observable (e.g., [2,3,4,6,15,16]). Some of the larger, and more areally extensive barranco systems have been forming over centuries or longer. Increased channel incision, along with increased sheet and rill erosion, is generally viewed as evidence of severe landscape degradation (sensu [17]). Yet, the creation or exacerbation of incision in barrancos is but one of many potential consequences of complex interactions between social and biophysical drivers of surface dynamics that have been shaping Mediterranean landscapes for millennia. The specific land-use histories of these coupled human and natural landscapes feedback into the earth-surface processes that shape them, in turn offering new constraints and opportunities for subsequent agropastoral and other land-uses [11]. It is therefore important to understand the potential long-term co-evolution between human land-use and barranco formation, however, the centuries-long time scales of these processes makes direct observation impossible. Proxy records of landscape change are widely and irregularly distributed in space and time, and often contain significant lacunae, allowing for multiple interpretations of the same evidence (e.g., [7,8]). Fortunately, advances in computational surface process modeling offer a way to investigate the complex, long-term interaction of anthropogenic and biophysical drivers of land-scape dynamics [18,19,20].

We describe the results of a series of modeling experiments, using a digital laboratory developed in the Mediterranean Landscape Dynamics (MedLanD) project, designed to explore the long-term consequences of small-holder agropastoral land-use for the evolution of barrancos in Mediterranean Spain. The MedLanD Modeling Laboratory (MML) is an open-source, integrated modeling environment that has the ability to couple spatially explicit (cellular automata) models of landscape evolution, agent-based and GIS-based models of human land-use, and regression or cellular models of vegetation and climate change to study the long-term dynamics of coupled human and natural landscapes [19,21,22,23,24,25,26]. In these experiments, we model the effects of increasing population, reducing fallowing intervals, and resource management strategies on barranco incision (Table 1). We situate these experiments in the real-world landscape of the Penaguila Valley in northern Alicante Province, Spain, which is the location of one of the earliest farming settlements (i.e., Neolithic) in the region [27,28]. Today, the valley is dissected by deeply entrenched barrancos containing incised sections that appear to postdate early Neolithic occupation of the valley.

Table 1. Table of modeling experiments conducted.

**Table 1.** Table of modeling experiments conducted.
Experiment Number	Number of People	Land Tenure Type
1	30	Satisfice
2	30	Maximize
3	60	Satisfice
4	60	Maximize
5	120	Satisfice
6	120	Maximize

2. MedLanD Modeling Laboratory (MML)

2.1. Surface Process Model

Many of the details of the surface process model component of the MML—r.landscape.evol—are described in Mitasova and colleagues [19]. (See Appendix Table A1 for a description of the input parameters of the module). In brief, r.landscape.evol is written in Python to run within the open-source GRASS GIS environment, where it can take advantage of fast computational hydrology tools, a parallelized map calculator, and special Python library. It uses a 3D implementation of the Unit Stream Power Erosion/Deposition (USPED) equation [29,30,31] to estimate transport capacity on hillslopes and rills (Equation (1)), and the Stream Power equation [19] to estimate transport capacity in channels (Equation (2)):

T_(hillslopes) = R K C A^m (sin β)ⁿ

(1)

T_(channels) = K_t n⁻¹ g_w h^m (tan β)ⁿ

(2)

where R, K, and C are the rain, soil erodibility, and land cover coefficients of the well-known RUSLE equation [32], A is the upslope accumulated area (per contour width), β is the local slope (in degrees), K_t is a coefficient of substrate erodibility in stream channels, n is Manning’s coefficient, g_w is the gravitational power of flowing water, h is the depth of flow, and m and n are empirically derived transport coefficients (both 1 for sheetwash on hillslopes, and 1.5 and 1.6, respectively, for flow in channels). The implementation of r.landscape.evol used here does not use soil creep or shear stress equations mentioned in Mitasova and colleagues [19], but these are alternative modes that exist in the module, and which may be implemented if desired.

The model estimates meters of erosion or deposition (ED) in each cell of a DEM on the basis of 2D divergence in transport capacity across topography (Equation (3)):

E D = \frac{d (T \cdot \cos (α))}{d x} + \frac{d (T \cdot \sin (α))}{d y}

(3)

where α is the topographical aspect (in degrees). A map of ED is converted to elevation changes by normalizing to sediment bulk density, and then added back to the DEM to create a new DEM to be the base layer for subsequent modeling. This process is iterated repeatedly to evolve the digital landscape. The model is transport capacity limited, but erosion is not unlimited. The model also requires a digital map of estimated bedrock elevations, and so tracks the amount of soil/regolith available for erosion. Erosion in excess of the available sediment in a given cell is not allowed. Although not implemented in the experiments presented here, bedrock erosion can be simulated by artificially deepening “soils” in places likely to experience significant erosion, and then “indurating” these soils with smaller values of K or K_t.

Sediment transport capacity can be altered by land cover, surface characteristics, or the amount of water available for runoff through the process equation coefficients [32,33] (Figure 1, see also Appendix Table A1). In the modeling experiments described here, K and R are kept constant, and take values empirically calculated for Mediterranean terra rossa soils [34] and mid-Holocene precipitation [22] (see Section 3.3, below). Human land-use activities, described below, can alter land cover. Land cover affects the calculation of erosion or deposition in two ways. On hillslopes the protective effect of vegetation on a plot of land is accounted for by C, affecting localized changes in transport capacity as estimated by the USPED equation. Additionally, the vegetation traps water, reducing runoff from a cell proportionally to the type of vegetation cover present (currently, runoff percentage is estimated from a linear regression of C vs. runoff water infiltration.). This changes the amount of water flowing though the downstream portions of the drainage, changing the estimated value of stream power in downstream reaches, and thereby affecting the calculation of erosion and deposition in those portions of the drainage. The input parameters of the surface process model are summarized in Appendix Table A1.

Finally, we parameterize the model to operate at a yearly interval to match better with the human agricultural cycle. This means that input parameters such as rainfall, storm frequency, vegetation growth, and land-use are all annualized, and that we do not explicitly model the occasional, very short term, extreme events that can have significant impacts on sediment transport (e.g., [14]). If these events were of significant effect in the formation of barrancos, then we may be underestimating the impact of human activity on barranco formation.

2.2. Human Land-Use Model

Although the MML contains a sophisticated Agent-Based land-use simulation engine, in this research we have chosen to use a more generalized model of Neolithic farming that simplifies land-use decisions with stochastic modeling techniques. This approach is more appropriate for modeling questions that do not require sophisticated agents (e.g., [22]), or that, as we do here, focus primarily on non-social aspects of socio-natural systems (but see [21,25,35] for modeling problems that do benefit from the ABM approach). Our simplified land-use model nevertheless encompasses a large range of human behaviors, however, and allows for several different land-use strategies to be modeled within a simple over-arching modeling framework that allows for faster simulation times.

Figure 1. A schematic representation of the effect of topographic and vegetation change on the calculation of erosion and deposition in the MML surface process model. Our modeling approach assumes that the sediment load in flowing water is normally near transport capacity (a spatially and temporally dynamic equilibrium influenced by factors like the quantity of water, slope, and land-cover). In a three-dimensional landscape, changes in one of these factors will alter transport capacity. Sediment will then be entrained or deposited until a local equilibrium between sediment load and transport capacity is reached. Hence, increases in slope (or, convexity) will tend to result in erosion; decreases in slope (or, concavity) will tend to result in deposition. Reduction of vegetation (shown as green bushes) will tend to result in increased erosion or decreased deposition, depending on the localized change in slope. (Length of blue arrows schematically represents transport capacity of overland flow).

The simplified land-use model amalgamates the decisions of individuals (or households) at the coarser social level of a village, which we here take to mean a group of individuals that cohabit in a community, and that generally communicate or coordinate subsistence and land-use decisions. The model is divided into three sequential components: (1) land-use planning; (2) enactment of a land-use plan; and (3) calculation of direct land-use impacts. The land-use plan for the village is updated on annual basis, and balances rudimentary information about mean productivity in the village’s catchments against a particular subsistence plan (ratio of agriculture to pastoralism), land choosing strategy (no land tenure with random allocation, satisficing subsistence needs, or maximizing agropastoral returns), and population, which are variables parameterized by the modeler at the start of a simulation run. Village labor availability and requirements, the general productivity and ecological characteristics of crops and herd animals, and the severity and spatial patterning of land-use impacts for herding or farming can all be parameterized to suit specific case-study conditions (see Appendix Table A2 for all land-use options of the model). Land-use occurs within predefined agricultural and pastoral catchments, which are created by the method described in [36], and which allows land above a slope threshold to be ignored (in our case, >10°). Each year, the village acquires new knowledge of the productivity of the landscape from the previous year’s returns. These data are added to a “memory bank” of variable length, and the subsistence plan for the upcoming year is created using the averaged values over a predetermined number of previous years. This information is subjected to a Gaussian random perturbation to prevent perfect knowledge of past events, which more realistically simulates the vagaries of human memory and human ability to estimate averages [37,38].

The village agent then uses this information to derive the number of farm plots and grazing patches that it believes are required to sustain the target population level (which does not change during the simulation). The land-use plan is enacted via spatially-explicit stochastic farming and grazing land-use choice procedures.

Three different procedures may be employed for farming land-use behavior, depending on the style of land-tenuring strategy desired for a particular model run. In the case of a non-tenuring farming strategy, new plots are randomly chosen from within the farming catchment every year. For a satisficing strategy, plots are randomly chosen the first year, and are never relinquished. If more plots are needed, they are randomly selected from the unused portion of the catchment as needed, and remain under cultivation for the remainder of the simulation. For a maximizing strategy, plots are also randomly chosen the first year, but may be dropped in subsequent years if their productivity falls below a predetermined proportion of the average yield. Randomly chosen new fields are then added as necessary in order to meet farming goals. Previously used-and-relinquished fields are available for reuse in the future. Farming returns are calculated per-plot by the average of three regression formulas, each calculated from empirical data about the effect of soil depth, soil fertility, and rainfall on wheat and barley yields in Mediterranean climate regimes, and calibrated to typical yields of ancient wheat and barley varieties [39,40,41,42,43,44,45,46,47,48].

Grazing occurs within a subset of the grazing catchment that is determined by the current number of grazing patches required and a least-cost routing algorithm. The number of grazing patches is calculated by the village agent each year based on the current village population, the parameterized herding ratio, the fodder requirements of herd animals, and the agent’s current knowledge of average grazing patch fodder yields. Optionally, grazing can be allowed on the fallowed portions of the farming catchment as well. Thus, non-denuded grazing patches near the village are always grazed, but farther patches may be added as necessary. Although all patches within the delimited subset of the grazing catchment will be grazed in a given year, some patches are grazed more intensively than others. The spatial patterning of this grazing intensity is determined by a Gaussian random function with tunable spatial autocorrelation that creates a patchy impact pattern. This more realistically simulates actual grazing patterns for site-tethered pastoralism [49,50,51]. Grazing returns are calculated from the known amount of digestible matter for given vegetation types, and the intensity of grazing [52,53,54,55,56,57,58,59]. People then gain calories from the milk and meat produced by the herds that were supported by the current grazing plan [60,61,62,63,64].

Once the amount and location of land-use activities are decided by the village agent for a given simulation-year, the model enacts the plan, assesses the amount of agricultural and pastoral returns that were gained, calculates the impacts of the land-use to vegetation and soil fertility, updates the landscape with the effect of these impacts, and then passes those new parameters to the surface process model. Impacts of farming a new plot include reduction of its vegetation cover to a value appropriate for cereals (i.e., dense grassland), and all farmed plots also experience reduction in soil fertility. Soil fertility is reduced according to a Gaussian random function with tunable spatial autocorrelation, with mean and standard deviation set at the start of the modeling run. Grazing reduces land-cover by the amount determined by the grazing intensity map discussed in the previous paragraph. Grazing also affects soil fertility in that manure from grazing animals is added to the soil at a rate commensurate to the intensity of grazing on that patch [65].

The map of land cover (converted to values of C) for a given year is passed as input into the r.landscape.evol surface process model as described above. The calculated erosion or deposition changes the surface topography, affecting the depth of soils and local slopes, which feeds back into the farming and grazing planning process for the next year.

3. Modeling Experiment Design

3.1. Creating the Digital Landscape

We focus the experiments for studying the long-term impacts of land-use practices on barranco formation in a digital representation of the real-world landscape of the Rio Penaguila Valley, near the town of Benilloba, Alicante (Figure 2). We use the early Neolithic archaeological site of Mas d’Is [27,66] as the location of a small-holder farmer/herder village. Today, this area is dissected by five large barrancos that probably did not exist, or were much less deeply incised, in the mid-Holocene. Hence we need to initiate our surface process models on a landscape that lacked these erosional features. To do so, we conducted geomorphological and geoarchaeological fieldwork in the project area, compiling a variety of temporal data for the different landforms to create a landscape chronology. These data included the date of the earliest archaeological material recovered during pedestrian surveys, stratigraphic information, and OSL dates where possible. We were able to delineate terraces and other areas that were likely present during the Neolithic period (“Terrace A” in Figure 3), as well as those areas that are more recent. In general, post-Neolithic surfaces were located in the bottoms of the barrancos and on low alluvial terraces in barranco-bottoms. We masked these more recent surfaces in a GIS to remove them from the digital landscape (i.e., DEM), and interpolated new terrain into the masked areas using the topographic information from the adjacent, older remnant landforms (e.g., “Terrace A”). The interpolation routine was tuned to also reduce the sharpness of slope curvatures (terrace edges), making slope breaks more natural. Any artificial internally-draining basins that were created in the valley floors by the interpolation routine were then filled so that the valley-bottoms were hydrologically continuous. This produced a landscape of broad U-shaped valleys, very different from the highly incised modern landscape (Figure 4).

Figure 2. Map of Southeast Spain, showing the location of the Rio Penaguila valley project area.

Figure 3. Reconstructed mid-Holocene topography in the Penaguila valley project area. The Neolithic site of Mas d’Is is shown, as are terraces near the site that are of Neolithic, or earlier, age (all labeled as “Terrace A”). The main watercourse passing by Mas d’Is is shown as a thick blue line, and the other portions of the extracted stream network are shown as thin blue lines. High elevations are shown in reds and brown, and lower elevations are shown in greens and yellows.

3.2. Modeling Past Climate

We retrodicted climate variables for the late Neolithic I period of southern Spain (7550-6450 cal. BP) with a regional downscaling method that localizes paleoclimate retrodictions from large-scale Global Climate Models (GCM) based on a regression relationship between the GCM output for the 30-year climate norm, and the actual observed climate information observed at a weather station [67,68,69]. We have used this method in several other research projects with the MML (e.g., [21,25,34,70]). We used the average climate values derived across the entire Neolithic I period, downscaled at a weather station in Benifallim, a town close to the Penaguila valley project area. The specific climate values that we used are shown in Appendix Table A1.

Figure 4. A comparison of the modern topography of a section of (a) the upper Rio Penaguila; and (b) the reconstructed Neolithic topography in the same area. The paleotopographic reconstruction resulted in a landscape of broad, U-shaped valleys, markedly different from the modern terrain.

3.3. Modeling Land-Use

We sited a small farming village at the location of the Mas d’Is archaeological site on the reconstructed mid-Holocene digital landscape, and performed a set of six experiments with the village population initialized at 30, 60, or 120 individuals using two different land-use strategies (Table 1). We constrained the area a village can use for cultivating cereal crops so that the villagers must reduce the fallow interval to produce sufficient food as the population grows. We also invoke two different strategies to guide agropastoral land use: satisficing and maximizing. With both strategies, a model begins by identifying all land that fits a set of basic criteria of suitability for cultivation and ovicaprine herding. In the experiments reported here, we parameterized the village based on ethnographic and archaeological data about small-scale, subsistence-focused, village-based Mediterranean agropastoralism. Neolithic subsistence in southern Europe was largely based on cereals [71]; therefore, we set 80% of a villages caloric needs to be met by consuming plants and 20% to be met by consuming animal products. Based on information about Neolithic agricultural practices [71,72,73], we parameterized 25% of the plant portion of the diet to be met by barley, and the remainder by wheat. Herds were parameterized as a 2:1 mix of goats and sheep, with the average herd animal requiring about 680 kg of fodder per year [36].

All together, this corresponds to a need of about six herd animals and about 230 kg of cereals per person per year. Herd animals were allowed to graze on fallowed portions of the agricultural catchment, as well as on the stubbles left after harvest of the farmed fields. Their manure contributed fertility gains of about 0.2% per year to the areas they grazed (including agricultural fields grazed for stubbles). Farm plots are modeled as discrete rectangular 20 m by 50 m regions, which are on par with the size of traditional, hand-worked historic peasant farm-plots in the southern Mediterranean region [74]. Farming reduced fertility at an average rate of 3% per year, and herd animals reduced land cover by an amount that would require about three years to regrow. Herds were parameterized as a 2:1 mix of goats and sheep, with the average herd animal requiring about 680 kg of fodder per year [36]. We kept the size and configuration of the farming and grazing catchments constant between all six experiments. This is not unreasonable, since real-world farming settlements are geographically constrained by the presence of other farmers in a region. Catchment modeling was conducted using a least-cost surfaces method described by Ullah [36]. For cereal cultivation, this included all land with slopes ≤ 10° within a 30 min walking distance of the village. Land suitable for herding was not limited by slope and extended to include all land within the Penaguila drainage (up to about an 8-h walk from the village).

Oversizing the farm catchment allows for the possibility of fallowing (particularly in maximizing strategies) on a time scale determined dynamically by the ratio of in-use fields to the catchment size. Initially, landscape cells are allocated to a village “agent” by randomly selecting enough suitable cells for cereal cultivation and enough suitable cells for animal pasturing to meet the caloric needs of the village (i.e., depending on the initial village population). The total number of fields under cultivation can change over time, depending on the productivity of previously used fields within the time span of the agent memory (five years). For example, with a village population of 120 and a 313 ha farming catchment (used in the experiments reported here), about 3.5% of the possible farming catchment, would be cultivated initially. But as farming returns from the initially-used fields declines, the same village will need to increase the farmed area to 6%–7% of the total to meet their needs. In the same way, the percentage of the grazing catchment under use can change over time. The same 120 people would use between 10% and 17% of the grazing catchment, depending upon their perception of the current grazing productivity of vegetation within the catchment.

With a “satisficing” strategy, a village will compare land in use with caloric needs each subsequent annual time-step of the model. If caloric needs are exceeded, some land (with the lowest returns from cereal cultivation and animal herding) will be fallowed. If caloric needs are not met, then additional suitable landscape cells will be selected randomly to be put into use. With a “maximizing” strategy, at each time-step, the village will evaluate the performance of all land with regard to its caloric returns. Underperforming land (that is, fields that produce less than 80% of the current average yield) will be relinquished to fallow and will be replaced with new land chosen randomly from among the suitable and currently unused landscape cells. Appendix Table A2 summarizes the range of values used to parameterized the social portions of our simulation experiment.

3.4. Sensitivity Testing and Experiment Repetition

The stochasticity embedded in modeling procedure, and the complex interactions between land-use and landscape change means that the results differ in each model run, even if the initial parameters are the same. To more accurately assess the impacts of agropastoral land-use on barranco erosion, it is necessary to repeat every experiment (i.e., model scenario with a particular set of initial parameters) multiple times, and calculate central tendencies and variance measures across multiple dimensions. While this does mean that the current research cannot focus on infrequent, “extreme” results (which may be important drivers of barranco evolution in some cases) the use of central tendencies allows us to understand the base-line (averaged) affects of human-land use on barranco formation under different land-use scenarios. The number of repetitions needed to adequately encompass the range of variation in land-use and landscape interactions can vary greatly with different model algorithms and model designs. Hence, we carried out a series of tests to assess model sensitivity to variation in initial parameters and estimate an optimum number of repetitions needed.

We carried out an initial run of the 60-person satisficing model, with 100 repetitions of 500 years each. We extracted several types of statistics from each model run (e.g., average amount of erosion/deposition per year, village population, proportion of different vegetation classes over time), and calculated the standard error of these statistics for all permutations of repetitions (1 through 100 repetitions). We repeated this process 100 times, shuffling and re-sampling the set of runs between each repeat using Monte Carlo methods. Example results of this sensitivity analysis show how the standard error for mean erosion/deposition in the stream channels changes as the number of run repetitions (n) increases (Figure 5). These tests indicated that error ceases to significantly improve after about 40 repetitions, and so we limited our other experiments to this number of repetitions. We also observed that the initial runs of 500 model years did not yield significantly different results from shorter runs. Additional experiments, therefore, ran for a total of 300 model years, significantly reducing processing-time and storage requirements. Running the model for 300 years still allowed for a period of about 50 years in which village land-use stabilized, followed by 250 years in which the land-use pattern of each experiment impacted erosion and deposition in the watershed. The combination of six experiments, repeated 40 times resulted in a total of 240 individual model runs of 300 years each (Table 1). These produced 40 time series of 300 maps for annual erosion/deposition, soil depth, land-cover, and topography for each experiment.

Because landscapes are also dynamic in the absence of any human presence, we conducted two additional experiments without any human land-use. In one experiment, we covered the landscape with “typical” Mediterranean woodland as would be present under mild climatic conditions with no farming or grazing; in the other one, we mantled it with sparse grass/herb cover as might be present under harsher climatic conditions. These provide points of comparison of human impacts with extremes of non-anthropogenic climate-induced landscape changes in this region.

For all experiments, we calculated total, cumulative maps of landscape change, as well as statistical maps of central tendencies and variance. The results of all repeated runs were statistically combined for each of the six experiments as maps and summary values, to allow comparison of the landscape consequences of different land-use configurations. We refer to these syntheses of landscape change in the discussions below. The results of these computational experiments in socio-ecological system dynamics are summarized below.

Figure 5. Results of sensitivity analysis about the effect of the number of experiment repetitions on the error surrounding mean estimates of erosion in the main Barranco near Mas d’Is. Red line is the mean of all recombinations at the specified number of repetition. The dark grey shaded region denotes 1 SD, and the lighter grey shaded region denotes 2 SD, based on the range of values from all 100 recombinations at each number of repetitions. The rate at which the standard error decreases with sample size (n) greatly tapers off after about 40 repetitions (indicated by the arrow).

4. Simulation Results

4.1. Erosion and Deposition in Barrancos

In order to use quantitative modeling to investigate the potential impacts of rural land-use on gullying and barranco incision it is necessary to have a surface process model that can accurately represent incision. While we have been developing the surface process model of the MML for nearly a decade, our focus was more on hillslope processes than incision [21,34,35]. Recently, we enhanced this model to better represent erosion and deposition in channels, combining well-known formulas to represent surface processes across the different landforms of small watersheds [19].

The results are encouraging. An oblique detailed view of the final topography for one run of one experiment, colored by the 300-year sum-total erosion/deposition for each map pixel, shows that the surface process model is capable of creating significant erosional dynamics in stream channels (Figure 6). Streams incise at the outer edge of channel bends, and at places where the channel gradient rapidly increases. Material is deposited where the gradient decreases and at the inside of channel bends. Flat areas are characterized by meanders and channel levies, and there are several remnants of abandoned channels throughout the stream course.

Figure 6. An oblique, three-dimensional view of a portion of the region shown in Figure 4. Warm colors (yellow through red) indicate cumulative erosion, and cool colors (green through blue) indicate cumulative deposition. The locations of the four cross-sections are indicated, as is the position of the Neolithic site of Mas D’is. Note that channel erosion occurs most intensively where narrowing of the channel bottom, coincides with a rapid convex change in slope in the channel gradient (e.g., points “c” and “d”). These are similar to the process of “head cuts” in real channels (e.g., [5]).

We sampled one of these digital landscapes through time at four cross sectional transects of the middle reaches of the major barranco that passes near Mas D’is, from its upstream to downstream reaches, indicated by the labeled horizontal lines in Figure 6. We recorded elevation cross-sections at these sampling transects at 50-year intervals in the 300 model-year time series at each of these transects (Figure 7). These cross-sections show the temporal dynamics of downcutting in different places along the drainage network. Barranco incision increases from upstream to downstream in the middle reaches of the barranco. This is to be expected as a result of increased water flow through the channel. Lateral incision is more prominent than absolute downcutting, but distinct channelization and incision did occur in all experiments. This is similar to empirical observations of the process of erosion in real barrancos [5]. Channel bar and levy development may offset some of the total erosion over time, however, particularly as the main channel of the barranco migrates within the valley-bottom. This indicates that incision and erosion of these barrancos may be a temporally complex phenomenon, with time-lag in sediment transport obscuring the observability of the total incision that is occurring in the system. This would mean that barranco downcutting would not have been immediately obvious to local farming peoples.

Figure 7. (a–d) Four cross-sections of the main Barranco: two profiles upstream from Mas d’Is, near Mas d’Is, and downstream from Mas d’Is. The red line shows the starting reconstructed mid-Holocene cross-section, and the green lines show the cross-section at 50-year intervals during one run of the maximizing large-population experiment. Incison and lateral erosion are apparent at all cross sections, although the absolute amount of incision is greater in the downstream cross sections. The locations of the cross-sections are shown in Figure 5 and Figure 8.

Over the total length of a barranco channel, however, erosion should decrease in the furthest downstream section, as local base-level is approached, to create a graded profile. We kept track of the longitudinal profile of the main Mas D’is barranco (indicated by the thick blue line in Figure 3) as it changed over the 300 model-years of the simulation experiments (Figure 8). The barranco profile elongates, and does become increasingly graded over time. The elongation is likely due to increased channel sinuosity as meanders form in the flats, and as lateral incision cuts into the valley walls. Overall, erosion tends to occur in punctuated stretches of the channel, producing a stepped longitudinal profile similar to head cuts in real-world gullies (Figure 5 and Figure 8).

Figure 8. The longitudinal profile of the main watercourse near Mas d’Is (see Figure 2). The red line indicates the profile as extracted from the starting reconstructed mid-Holocene topography. The green lines show the profile in 50-year increments during one run of the maximizing large-population experiment. The overall length of the profile increases over time due to changes in channel morphology (such as the growth of meanders and lateral incision). The profile also becomes more graded over time, although erosion seems to occur mainly in punctuated stretches of the channel bottom.

The range of values for average channel incision after 300 model-years varies from 0.062 to 0.083 m (or, 0.21 to 0.28 mm/yr), with a mean value for average channel incision of 0.07 m (0.23 mm/yr) for all experiments combined. These rates are similar to those observed for recent erosion in barrancos [3,15,16]. In each of the experiments, there were locations in the channel that were eroded to depths of 5 to 6 m at the end of the simulation period, again similar to rates of head cutting in modern barrancos [3,15,16]. To be clear, although we initialized the MML surface process model with realistic values, we have not yet attempted to fine-tune the internal model dynamics to match those of real-world landscapes. In part, this is a function of the lack of detailed, quantitative data on barranco erosion rates at century to millennial scales. What we have tested here is whether our modeling environment can simulate real-world erosional processes, even if the simulations of specific quantities of sediment moved may not match those of Mediterranean barrancos in the Holocene. In fact, the tests summarized here indicate that the surface process component of the MML creates realistic patterns of channel incision in both the longitudinal and lateral dimensions of barrancos. This improves our confidence in the results when we use this environment to examine interactions between human land-use and landscape change. We review these socio-ecological interactions below.

4.2. Effects of Land-Use on Barranco Incision

As noted above, we evaluated the consequences for landscape change of six combinations of population size and decision strategies (Table 1). We parsed the landscape into three socio-ecological meaningful units: the Mas D’is barranco valley-bottom (Figure 9, transparent outlined area), the portions of the landscape outside of this barranco (mostly hillslopes and smaller channels), and the modeled Mas D’is farming catchment (Figure 9, grey outlined area). Partitioning the landscape in this way helps us to better understand the connection between human land use and the spatial patterning of erosive impacts. We do this by separately calculating the cumulative amount of sediment eroded and deposited in each of these three regions after 300 simulated years in each experiment.

Figure 9. Map of cumulative erosion and deposition after 300 simulated years for one run of the maximizing large population scenario (see Table 1). The farming catchment is shown as a shaded region surrounding the site of Mas d’Is (indicated by the red star). The outline of the main Barranco that passes near Mas d’Is is also shown, as are the locations of the four cross-sections shown in Figure 6. Warm colors (yellow through red) indicate cumulative erosion, and cool colors (green through blue) indicate cumulative deposition. Scale is in meters.

At the scale of the watershed (Figure 10), 300 years of human farming and animal herding produces erosion rates on hillslopes and small channels that slightly exceed those for the natural wooded landscape (transparent yellow line), but are considerably less than erosion in natural sparse, open vegetation (solid yellow line). The difference between these two control runs simulates large-scale shifts in vegetation cover that would accompany climate change, and provide interesting benchmarks to compare to the simulations that included agropastoral activity. Sediment deposited on hillslopes in human-occupied landscapes is roughly the same as in an unoccupied wooded landscape. Overall, small-holder agropastoralism has limited effects on hillslope dynamics, at least initially. Unsurprisingly, the simulated community with the largest population and following a returns-maximizing strategy had the strongest impact, but only by a comparatively small amount. Population makes a slightly larger difference than satisficing vs. maximizing strategy.

Figure 10. Cumulative amount of sediment eroded or deposited (tons/ha) on hillslopes for each of the six experiments and two control runs during the 300-year modeling interval. Values are averaged across the 40 repetitions of each experiment. The cumulative tons/ha of eroded sediment are shown below the center lines of each graph, and cumulative tons/ha of deposited sediment are above the center lines for the 300 year modeling interval. We show erosion and deposition separately because sediment removal and accumulations often happen in different landscape contexts, with different consequences for agropastoral and other land-use.

Within the farming catchment (Figure 11), anthropogenic erosion from all scenarios exceeds that of natural vegetation change, while anthropogenically-driven deposition falls between the wooded and open landscapes. Interestingly, a maximizing strategy and medium sized population (solid blue line) produces considerably more erosion than other agropastoral land-use scenarios, including those with higher populations. In other words, small-holder agropastoral land-use has significant impact at local scales, but is less important than large-scale vegetation change (e.g., due to climate change). Additionally, while the degree of impact seems to scale with the number of people and their subsistence activities at the watershed scale, at the local level, the strategy used to make land-use decisions can exceed the landscape impacts of increased population density alone.

Figure 11. Cumulative amount of sediment eroded or deposited (tons/ha) in the delineated Mas d’Is farming catchment (see Figure 4) for each of the six experiments and two control runs during the 300-year modeling interval. Values are averaged across the 40 repetitions of each experiment. The cumulative tons/ha of eroded sediment are shown below the center lines of each graph, and cumulative tons/ha of deposited sediment are above the center lines for the 300 year modeling interval. We show erosion and deposition separately because sediment removal and accumulations often happen in different landscape contexts, with different consequences for agropastoral and other land-use.

The natural and human effects on barranco incision are much different from those on the rest of the watershed. Whereas a large-scale shift from natural woodland to open vegetation on the hillslopes and small channels increases both erosion and deposition by roughly equal amounts within the watershed (Figure 10), within the barranco valley, non-anthropogenic vegetation change has disproportionate effects on erosion and deposition. That is, at the temporal scale of 300 model years, the effect of large-scale denudation is to increase the overall sedimentary dynamics of hillslopes, but not in the main channel. Erosion does increase slightly in the main barranco valley (at a level more or less equivalent to that in the rest of the watershed) but deposition increases by several orders of magnitude. The major effects of climate-induced large-scale devegetation, then, would be hillslope erosion and valley alluviation, agreeing with observations of field geomorphologists (e.g., [75]).

Anthropogenic impacts on the main barranco differ significantly from non-human vegetation changes (Figure 12). Starting with a wooded landscape, both erosion and deposition are greatly increased by small-holder agropastoral land-use of all kinds modeled here. Anthropogenic deposition rates exceed those of non-human open woodland, and anthropogenic erosion rates in barrancos exceed by order of magnitude those driven by non-human vegetation change, as well as all erosion rates in the rest of the watershed. So, human activity contributes little in the way of sediment to these watercourses, but greatly increases the dynamics of sediments already in the floors and banks of barrancos. This effect seems more strongly affected by the number of people engaged in agropastoral activities than the strategy they use, at least for erosion. However, the overall rate of erosion slows over time in the anthropogenic cases, while the rate of deposition does not. This may indicate that incision in barranco systems begins to occur relatively quickly after initial anthropogenic land cover changes, but, without further changes in land-use or climate, may eventually achieve a metastable equilibrium (balance of erosion and deposition).

Figure 12. Cumulative amount of sediment eroded or deposited (tons/ha) in the Mas d’Is Barranco for each of the six experiments and two control runs during the 300-year modeling interval. Values are averaged across the 40 repetitions of each experiment. The cumulative tons/ha of eroded sediment are shown below the center lines of each graph, and cumulative tons/ha of deposited sediment are above the center lines for the 300 year modeling interval. We show erosion and deposition separately because sediment removal and accumulations often happen in different landscape contexts, with different consequences for agropastoral and other land-use.

5. Conclusions

The MedLand Modeling Laboratory provides a useful tool to examine surface dynamics and the impacts of rural land-use on landscapes [19,21,34]. The inclusion of separate hillslope- and stream-process models allow it to achieve reasonable results on slopes, while also more realistically simulating erosion and deposition in intermittent watercourses like barrancos that dissect landscapes of the western Mediterranean. The results presented here show that the MML can simulate simultaneous vertical and lateral erosion, generate head-cuts, and produce graded profiles over time (Figure 6, Figure 7 and Figure 8). This allows us to investigate the impacts of the spread of agropastoral communities and associated land-use practices in the western Mediterranean, and clarifies the potential role of humans in barranco formation. Whereas our simulation experiments show that large-scale (climate-driven) vegetation change does greatly increase incision in these watercourses, they also show that smaller-scale, localized vegetation changes brought about by human land-use can have similar effects on barranco erosion rates, without climate change (Figure 10, Figure 11 and Figure 12).

McClure and colleagues [76,77] suggest on the basis of archaeological evidence that the earliest Neolithic agropastoralists in Mediterranean Spain preferentially cultivated the alluvial soils of valley bottoms, but subsequently shifted to upland localities. They hypothesize that a combination of expanding agropastoral land-use and changing climate made these valley-bottom locales vulnerable to severe erosion, forcing a reorganization of settlement and agricultural practices that appear archaeologically as the Late Neolithic (Neolithic II-b locally). Using an earlier version of the MML, Barton and colleagues [34,35] show how modest increases in the size of small-holder agropastoral communities in northern Jordan, engaged in wadi-bottom farming, could have triggered a phase shift in erosion/deposition ratios, leading to declining agricultural productivity.

The modeling results presented here appear to confirm these earlier studies and offer a more nuanced view of the relationships between human activity, climate change, and the location and intensity of landscape change. Climate-driven decrease in vegetation cover can increase the surface dynamics and sediment transport rates on hillslopes and within watersheds (Figure 10), with increased alluviation in intermittent watercourses like barrancos and wadis (Figure 12). While small-holder agropastoralism can have significant local effects on hillslopes within the catchment actively cultivated (Figure 11), its impacts on the broader landscape are minimal (Figure 10)—at least until farms and pasturage cover significant proportions of a region. However, even very limited settlement and agropastoral land-use appears to have impacts on the local and downstream dynamics of intermittent watercourses that are orders of magnitude greater than non-anthropogenic vegetation change at regional scales (Figure 12). The greatly increased sedimentary dynamics in the farming catchment erodes farmed fields in some locales and buries them in others, with the potential for significant loss of productivity and greatly increased risk for subsistence farmer/herders (both in terms of uncertainty and potential for loss). The easily-cultivated and fertile soils of these watercourses, which are so common in Mediterranean landscapes, would have made them especially attractive to the earliest farmers of the region, who relied on hoes and digging sticks instead of plows. Cultivating these valley-bottoms, however, could lead to their increasing unsuitability for agriculture within a few generations—potentially with correlated impacts spanning areas much larger than just that of the farming communities or their catchments of active use. The increased downstream incision that results from valley-bottom farming would decrease the number of viable fields in the alluvial zone, limiting future options for valley-bottom cultivation when the fertility of the local catchment is depleted, or when local agricultural soils have themselves been severely eroded. McClure and colleagues suggest that the loss of these initially highly productive patches of the Mediterranean landscape created strong socioeconomic incentives for significant reorganization of both subsistence practices and society toward greater complexity [76,77].

The kinds of processes that we have presented here also exemplify the kinds of non-linear dynamics that typify human-environmental interactions in complex, coupled socio-natural systems. The difficulty of predicting these non-linear dynamics through normal reductionist science or empirical statistical analyses underscores the importance of the kinds of modeling illustrated in this Special Issue of Land. Computational modeling allows investigation of feedbacks between the human and natural components of earth-systems processes with reproducible and controlled experiments that have tangible implications for the way we understand the empirical record of landscape change. We wish to stress that computational and empirical research techniques complement, rather than compete with, each other, and integrating the two approaches is essential for a holistic approach to understanding earth-surface change. For example, the simulations discussed in this paper were designed to extend information gained from empirical observation of human land-use and barranco erosion over short periods to gain insight into socio-ecological processes that play out over centuries. Ultimately, we learn not only that relatively moderate changes to vegetation by small populations of agropastoral villagers can greatly affect the capacity for incision in barrancos, but that very deep barranco incision may require continual changes to land-use (or climate) over long periods of time. It may also be that self-amplifying feedback processes between land-use and erosion of easily cultivated land in barranco bottoms could have driven the kinds of increases in social complexity and land-management that would have continually reinforced erosive regimes, producing the deeply incised landscapes that exist throughout southern Europe today. These insights were made possible by the combination of traditional and simulation research approaches.

Acknowledgments

This work was supported by a grant from the US National Science Foundation, Coupled Natural and Human Systems Program, Grant #DEB-1313727.

Author Contributions

C. Michael Barton directs the Mediterranean Landscape Dynamics project, planned and directed this research component of the larger project, and wrote the majority of the paper. Isaac Ullah carried out the modeling experiments and wrote a significant portion of the paper. Arjun Heimsath advised Barton and Ullah on geomorphic dyamics modeled in this research and wrote sections of the paper.

Appendix

Table A1. Table of landscape evolution parameters used in the surface process modeling component of the MML.

**Table A1.** Table of landscape evolution parameters used in the surface process modeling component of the MML.
Landscape Evolution Parameter	Typical Range of Values	Units	Value(s) Used in This Paper	Explanation
USPED K-factor	0.01–0.75	((T·ha·hr)/(ha·MJ·mm))	0.42	Soil erodibility index from RUSLE.
USPED R-factor	0–50	((MJ·mm)/(ha·hr·yr))	4.54	Rainfall factor from RUSLE.
USPED C-factor	0.005–0.5	unitless	0.005–0.5	Vegetation cover factor from RUSLE.
Kt	0.001–0.000001	unitless	0.0001	Stream transport efficiency variable (erodibility of stream substrate).
Sediment load exponent	1.5,2.5	unitless	1.5	Stream transport type variable (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport).
Manning's N	0.01–0.16	s/m^1/3	0.05	Average value of Manning's surface roughness coefficient value for channelized flow in the drainage.
Flow speed	0-3	m/s	1.4	Average velocity of flowing water in the drainage.
Soil density	0-3	T/m³	1.2184	Soil density map or constant for conversion from mass to volume.
Transition point	10–500	number of raster cells	100	Flow accumulation breakpoint value for shift from hillslopes to stream flow.
Per-storm precipitation totals	0–10,000	mm	20.61	Precipitation totals for the average storm.
Number of storms	0–300	storms/yr	25	Average number of storms per year.
Storm length	0–72	hr	24	Length of the average storm.

Table A2. Table of agropastoral parameters used in the stochastic agropastoral simulation component of the MML.

**Table A2.** Table of agropastoral parameters used in the stochastic agropastoral simulation component of the MML.
Village Subsistence Characteristic	Typical Range of Values	Units	Value(s) Used in This Paper	Explanation
General Village Characteristics
Number of people in the village	0–1000	number of individuals	30, 60, 120	The “target” number of people to be fed every year. Stays constant throughout the simulation.
Length of village “memory”	0–59	yr	5	Length of the “memory” of the agent in years. The agent will use the mean surplus/deficit information from this many of the most recent previous years when making a subsistence plan for the current year.
Amount of agricultural labor available	0–365	person-days	300	The amount of agricultural labor an average person of the village can do in a year.
Required amount of cereals	300–500	kg	370	Amount of cereals that would be required per person per year if cereals were the only food item being consumed.
Required number of animals	40–100	number of individuals	60	Number of herd animals that would be needed per person per year if pastoral products were the only food item being consumed.
Required amount of animal fodder	500–1000	kg	680	Amount of fodder required per herd animal per year.
Agropastoral Ratio	0–1	unitless ratio	0.2	Actual ratio of agricultural to pastoral foods in the diet, where 0 = 100% agricultural and 1 = 100% pastoral.
Village Farming Characteristics
Agricultural mix	0–1	unitless ratio	0.25	The wheat/barley ratio (e.g., 0.0 for all wheat, 1.0 for all barley, 0.5 for an equal mix).
Field dimensions	5–100	m	20, 50	North-South and East-West dimensions of agricultural fields.
Labor per field	5–100	person-days	50	Number of person-days required to till, sow, weed, and harvest one farm field in a year.
Field landcover value	0–50	succession stage	5	The landcover value for farmed fields (corresponds to an appropriate value from the landcover regrowth scheme).
Farming impact	0–10 (0–5)	% of maximum fertility	3 (2)	The mean and standard deviation of the amount to which farming a patch decreases its fertility (in percentage points of maximum fertility). Fertility impact values of individual farm plots is randomly chosen from a gaussian distribution that has this mean and standard deviation.
Maximum wheat	3000–4000	kg/ha	3500	Maximum amount of wheat that can be grown.
Maximum barley	2000–3000	kg/ha	2500	Maximum amount of barley that can be grown.
Satisficing farming strategy	Y/N	boolean	Both Y and N	Land is never dropped, only added if needed.
Maximizing Farming strategy	Y/N	boolean	Both Y and N	Land is dropped if below a previously defined threshold in productivity.
Productivity threshold	0–1	unitless ratio	0.2	Threshold for dropping land out of tenure with a maximizing strategy, interpreted as a percentage below the yearly average yield of all farm cells.
Fertility regain rate	0–100 (0–100)	% of maximum fertility	2 (0.5)	The mean and standard deviation of the natural fertility recovery rate (percentage by which soil fertility increases per year if not farmed). Fertility recovery values of individual landscape patches will be randomly chosen from a gaussian distribution that has this mean and standard deviation.
Village Grazing Characteristics
Minimum grazability	0–50	succession stage	2	Minimum amount of vegetation on a cell for it to be considered grazable by ovicaprines (corresponds to an appropriate value from the landcover regrowth scheme).
Grazing spatiality coefficient	0–200	m	50	Spatial dependency of the grazing pattern in map units. This value determines how “clumped” grazing patches will be. A value close to 0 will produce a perfectly randomized grazing pattern with patch size equal to raster cell resolution, and larger values will produce increasingly clumped grazing patterns, with the size of the patches corresponding to the value given.
Grazing patchiness coefficient	0–1	unitless	1	Coefficient that, along with the spatiality coefficient, determines the patchiness of the grazing pattern. Value must be non-zero, and usually will be ≤1.0. Values close to 0 will create a patchy grazing pattern, values close to 1 will create a "smooth" grazing pattern. Actual grazing patches will be sized to the resolution of the input landcover map.
Maximum grazing impact	0–50	succession stage	3	Maximum impact of grazing in units of “landcover succession” per annual grazing event. Grazing impact values of individual patches will be chosen from a gaussian distribution between 1 and this maximum value (i.e., most values will be between 1 and this value). Value must be ≥1.
Manuring rate	0–100	% of maximum fertility	0.2	Base rate that animal dung contributes to fertility increase on a grazed patch in units of percentage of maximum fertility regained per increment of grazing impact. Actual fertility regain values are thus calculated as “manuring rate x grazing impact”, so this variable interacts with the grazing impact settings.
Avoid grazing in agricultural catchment	Y/N	boolean	N	If turned on, ovicaprines will not graze in unused portions of the agricultural catchment (i.e., do not graze on "fallowed" fields, and thus no “manuring” of those fields will occur).
Avoid grazing on field stubbles	Y/N	boolean	N	If turned on, ovicaprines will not do any “stubble grazing” on harvested fields (and thus no “manuring” of fields).

Conflicts of Interest

The authors declare no conflict of interest.

References

Nogueras, P.; Burjachs, F.; Gallart, F.; Puigdefàbregas, J. Recent gully erosion in the El Cautivo badlands (Tabernas, SE Spain). Catena 2000, 40, 203–215. [Google Scholar] [CrossRef]
Casalı́, J.; López, J.J.; Giráldez, J.V. Ephemeral gully erosion in southern Navarra (Spain). Catena 1999, 36, 65–84. [Google Scholar] [CrossRef]
Valcárcel, M.; Taboada, M.T.; Paz, A.; Dafonte, J. Ephemeral gully erosion in northwestern Spain. Catena 2003, 50, 199–216. [Google Scholar] [CrossRef]
Gutiérrez, Á.G.; Schnabel, S.; Contador, F.L. Gully erosion, land use and topographical thresholds during the last 60 years in a small rangeland catchment in SW Spain. Land Degrad. Dev. 2009, 20, 535–550. [Google Scholar] [CrossRef]
Marzolff, I.; Ries, J.B.; Poesen, J. Short-term versus medium-term monitoring for detecting gully-erosion variability in a Mediterranean environment. Earth Surf. Process. Landf. 2011, 36, 1604–1623. [Google Scholar] [CrossRef]
Hooke, J.M. Human impacts on fluvial systems in the Mediterranean region. Geomorphology 2006, 79, 311–335. [Google Scholar] [CrossRef]
Van Andel, T.H.; Zanagger, E. Landscape stability and destabilization in the prehistory of Greece. In Man’s Role in the Shaping of the Eastern Mediterranean Landscape; Bottema, S., Entjes-Nieborg, G., van Zeist, W., Eds.; A.A. Balkema: Rotterdam, The Netherlands, 1990; pp. 139–157. [Google Scholar]
Bintliff, J. Time, process and catastrophism in the study of Mediterranean alluvial history: A review. World Archaeol. 2002, 33, 417–435. [Google Scholar] [CrossRef]
Redman, C.L.; Fish, P.R.; James, S.R.; Rogers, J.D. The Archaeology of Global Change: The Impact of Humans on Their Environment; Smithsonian Books: Washington, DC, USA, 2004. [Google Scholar]
Perevolotsky, A.; Seligman, N.G. Role of grazing in Mediterranean rangeland ecosystems. BioScience 1998, 48, 1007–1017. [Google Scholar] [CrossRef]
Barton, C.M.; Bernabeu Auban, J.; Garcia Puchol, O.; Schmich, S.; Molina Balaguer, L. Long-term socioecology and contingent landscapes. J. Archaeol. Method Theory 2004, 11, 253–295. [Google Scholar] [CrossRef]
Fumanal Garcia, M.P. Dinámica sedimentaria Holocena en valles de cabecera del País Valenciano. Cuatern. Geomorfol. 1990, 4, 93–106. [Google Scholar]
Hill, J.B. Land use and an archaeological perspective on socio-natural studies in the Wadi Al-Hasa, West-Central Jordan. Am. Antiq. 2004, 69, 389–412. [Google Scholar] [CrossRef]
García-Ruiz, J.M. The effects of land uses on soil erosion in Spain: A review. Catena 2010, 81, 1–11. [Google Scholar] [CrossRef]
Martı́nez-Casasnovas, J.A. A spatial information technology approach for the mapping and quantification of gully erosion. Catena 2003, 50, 293–308. [Google Scholar] [CrossRef]
Vandekerckhove, L.; Poesen, J.; Govers, G. Medium-term gully headcut retreat rates in Southeast Spain determined from aerial photographs and ground measurements. Catena 2003, 50, 329–352. [Google Scholar] [CrossRef]
Hill, J.B. What difference does environmental degradation make? In The Archaeology of Environmental Change; Fisher, C.T., Hill, J.B., Feinman, G.M., Eds.; The University of Arizona Press: Tucson, AZ, USA, 2009; pp. 160–173. [Google Scholar]
Clevis, Q.; Tucker, G.E.; Lock, G.; Lancaster, S.T.; Gasparini, N.; Desitter, A.; Bras, R.L. Geoarchaeological simulation of meandering river deposits and settlement distributions: A three-dimensional approach. Geoarchaeology 2006, 21, 843–874. [Google Scholar] [CrossRef]
Mitasova, H.; Barton, C.M.; Ullah, I.I.T.; Hofierka, J.; Harmon, R.S. GIS-based soil erosion modeling. In Treatise in Geomorphology: Vol. 3 Remote Sensing and GI Science in Geomorphology; Shroder, J., Bishop, M., Eds.; Academic Press: San Diego, CA, USA, 2013; pp. 228–258. [Google Scholar]
Wainwright, J. Can modelling enable us to understand the rôle of humans in landscape evolution? Geoforum 2008, 39, 659–674. [Google Scholar] [CrossRef]
Barton, C.M.; Ullah, I.I.T.; Bergin, S.M.; Mitasova, H.; Sarjoughian, H. Looking for the future in the past: Long-term change in socioecological systems. Ecol. Model. 2012, 241, 42–53. [Google Scholar] [CrossRef]
Barton, C.M.; Ullah, I.I.T.; Mitasova, H. Computational modeling and Neolithic socioecological dynamics: A case study from Southwest Asia. Am. Antiq. 2010, 75, 364–386. [Google Scholar] [CrossRef]
Mayer, G.R.; Sarjoughian, H.S. Composable cellular automata. Simulation 2009, 85, 735–749. [Google Scholar] [CrossRef]
Mayer, G.R.; Sarjoughian, H.S.; Allen, E.K.; Falconer, S.E.; Barton, C.M. Simulation modeling for human community and agricultural landuse. In Agent-Directed Simulation, Proceedings of the Agent-Directed Simulation Multi-Conference, Huntsville, AL, USA, 2–6 April 2005; Society for Computer Simulation International: San Diego, CA, USA, 2006; pp. 65–72. [Google Scholar]
Ullah, I.I. T.; Bergin, S. Modeling the consequences of village site location: Least cost path modeling in a coupled GIS and agent-based model of village agropastoralism in eastern Spain. In Least Cost Analysis of Social Landscapes: Archaeological Case Studies; White, D.A., Surface-Evans, S.L., Eds.; University of Utah Press: Salt Lake City, UT, USA, 2012; pp. 155–173. [Google Scholar]
Barton, C.M.; Ullah, I.I.; Mayer, G.R.; Bergin, S.M.; Sarjoughian, H.S.; Mitasova, H. MedLanD Modeling Laboratory v.1. CoMSES Computational Model Library. Available online: https://www.openabm.org/model/4609/version/1 (accessed on 23 April 2015).
Bernabeu Auban, J.; Orozco Köhler, T.; Diez Castillo, A.; Gomez Puche, M. Mas d’Is (Penàguila, Alicante): Aldeas y recintos monumentales del Neolítico Antiguo en el Valle del Serpis. Trab. Prehist. 2003, 60, 39–59. [Google Scholar]
Bernabeu, J.; García Puchol, O.; Pardo, S.; Barton, M.; McClure, S.B. AEA 2012 Conference reading: Socioecological dynamics at the time of Neolithic transition in Iberia. Environ. Archaeol. 2014, 19, 214–225. [Google Scholar] [CrossRef]
Mitas, L.; Mitasova, H. Distributed soil erosion simulation for effective erosion prevention. Water Resour. Res. 1998, 34, 505–516. [Google Scholar] [CrossRef]
Mitasova, H.; Hofierka, J.; Zlocha, M.; Iverson, R. Modeling topographic potential for erosion and deposition using GIS. Int J. Geogr. Inf. Syst. 1996, 10, 629–641. [Google Scholar] [CrossRef]
Moore, I.D.; Burch, G.J. Physical basis of the length-slope factor in the Universal Soil Loss Equation. Soil Sci. Soc. Am. J. 1986, 50, 1294–1298. [Google Scholar] [CrossRef]
Renard, K.G.; Foster, G.R.; Weesies, G.A.; Porter, J.P. RUSLE: Revised Universal Soil Loss Equation. J. Soil Water Conserv. 1991, 46, 30–33. [Google Scholar]
Renard, K.G.; Foster, G.R.; Weesies, G.A.; McCool, D.K.; Yoder, D.C. Predicting soil erosion by water: A guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). In Agriculture Handbook; US Department of Agriculture: Washington, DC, USA, 1997; Volume 703, pp. 1–251. [Google Scholar]
Onori, F.; de Bonis, P.; Grauso, S. Soil erosion prediction at the basin scale using the Revised Universal Soil Loss Equation (RUSLE) in a catchment of Sicily (southern Italy). Environ. Geol. 2006, 50, 1129–1140. [Google Scholar] [CrossRef]
Barton, C.M.; Ullah, I.I.; Bergin, S. Land use, water and Mediterranean landscapes: Modelling long-term dynamics of complex socio-ecological systems. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 2010, 368, 5275–5297. [Google Scholar] [CrossRef] [PubMed]
Ullah, I.I.T. A GIS method for assessing the zone of human-environmental impact around archaeological sites: A test case from the Late Neolithic of Wadi Ziqlâb, Jordan. J. Archaeol. Sci. 2011, 38, 623–632. [Google Scholar] [CrossRef]
Koriat, A.; Goldsmith, M.; Pansky, A. Toward a psychology of memory accuracy. Annu. Rev. Psychol. 2000, 51, 481–537. [Google Scholar] [CrossRef] [PubMed]
Schacter, D.L. The seven sins of memory: Insights from psychology and cognitive neuroscience. Am. Psychol. 1999, 54, 182–203. [Google Scholar] [CrossRef] [PubMed]
Quiroga, A.; Funaro, D.; Noellemeyer, E.; Peinemann, N. Barley yield response to soil organic matter and texture in the Pampas of Argentina. Soil Tillage Res. 2006, 90, 63–68. [Google Scholar] [CrossRef]
Araus, J.L.; Febrero, A.; Buxó, R.; Camalich, M.D.; Martin, D.; Molina, F.; Rodriguez-Ariza, M.; Romagosa, I. Changes in carbon isotope discrimination in grain cereals from different regions of the western Mediterranean Basin during the past seven millennia. Palaeoenvironmental evidence of a differential change in aridity during the late Holocene. Glob. Change Biol. 1997, 3, 107–118. [Google Scholar] [CrossRef]
Araus, J.L.; Amaro, T.; Zuhair, Y.; Nachit, M.M. Effect of leaf structure and water status on carbon isotope discrimination in field-grown durum wheat. Plant Cell Environ. 1997, 20, 1484–1494. [Google Scholar] [CrossRef]
Barzegar, A.R.; Yousefi, A.; Daryashenas, A. The effect of addition of different amounts and types of organic materials on soil physical properties and yield of wheat. Plant Soil 2002, 247, 295–301. [Google Scholar] [CrossRef]
Carter, D.L.; Berg, R.D.; Sanders, B.J. The effect of furrow irrigation erosion on crop productivity. Soil Sci. Soc. Am. J. 1985, 49, 207–211. [Google Scholar] [CrossRef]
Pswarayi, A.; van Eeuwijk, F.A.; Ceccarelli, S.; Grando, S.; Comadran, J.; Russell, J.R.; Francia, E.; Pecchioni, N.; Li Destri, O.; Akar, T.; et al. Barley adaptation and improvement in the Mediterranean basin. Plant Breed. 2008, 127, 554–560. [Google Scholar] [CrossRef]
Sadras, V.O.; Calvino, P.A. Quantification of grain yield response to soil depth in soybean, maize, sunflower, and wheat. Agron. J. 2001, 93, 577. [Google Scholar] [CrossRef]
Araus, J.L.; Slafer, G.A.; Romagosa, I.; Molist, M. FOCUS: Estimated wheat yields during the emergence of agriculture based on the carbon isotope discrimination of grains: Evidence from a 10th millennium BP site on the Euphrates. J. Archaeol. Sci. 2001, 28, 341–350. [Google Scholar] [CrossRef]
Wong, M.T.F.; Asseng, S. Yield and environmental benefits of ameliorating subsoil constraints under variable rainfall in a Mediterranean environment. Plant Soil 2007, 297, 29–42. [Google Scholar] [CrossRef]
Slafer, G.A.; Romagosa, I.; Araus, J.L. Durum wheat and barley yields in antiquity estimated from 13C discrimination of archaeological grains: A case study from the western Mediterranean Basin. Funct. Plant Biol. 1999, 26, 345–352. [Google Scholar] [CrossRef]
Adler, P.; Raff, D.; Lauenroth, W. The effect of grazing on the spatial heterogeneity of vegetation. Oecologia 2001, 128, 465–479. [Google Scholar] [CrossRef]
Alados, C.L.; ElAich, A.; Papanastasis, V.P.; Ozbek, H.; Navarro, T.; Freitas, H.; Vrahnakis, M.; Larrosi, D.; Cabezudo, B. Change in plant spatial patterns and diversity along the successional gradient of Mediterranean grazing ecosystems. Ecol. Model. 2004, 180, 523–535. [Google Scholar] [CrossRef]
Alados, C.L.; Pueyo, Y.; Barrantes, O.; Escós, J.; Giner, L.; Robles, A.B. Variations in landscape patterns and vegetation cover between 1957 and 1994 in a semiarid Mediterranean ecosystem. Landsc. Ecol. 2004, 19, 543–559. [Google Scholar] [CrossRef]
Al-Jaloudy, M.A. Country Pasture/Forage Resource Profiles: JORDAN; Food and Agriculture Organization of the United Nations: Rome, Italy, 2006. [Google Scholar]
Gulelat, W. Household Herd Size among Pastoralists in Relation to Overstocking and Rangeland Degradation (Sesfontein, Namibia). Master’s Thesis, International Institute for Geo-Information Science and Earth Observations, Enschede, The Netherlands, 2002. [Google Scholar]
Lubbering, J.M.; Stuth, J.W.; Mungall, E.C.; Sheffield, W.J. An approach for strategic planning of stocking rates for exotic and native ungulates. Appl. Anim. Behav. Sci. 1991, 29, 483–488. [Google Scholar] [CrossRef]
Nablusi, H.; Ali, J.M.; Abu Nahleh, J. Sheep and Goat Management Systems in Jordan: Traditional and Feedlot—A Case Study; Task Force Documents; Amman, Jordan, 1993. [Google Scholar]
Ngwa, A.T.; Pone, D.K.; Mafeni, J.M. Feed selection and dietary preferences of forage by small ruminants grazing natural pastures in the Sahelian zone of Cameroon. Anim. Feed Sci. Technol. 2000, 88, 253–266. [Google Scholar] [CrossRef]
Hocking, D.; Mattick, A. Dynamic Carrying Capacity Analysis as Tool for Conceptualising and Planning Range Management Improvements, with a Case Study from India; Overseas Development Institute, Pastoral Development Network: London, UK, 1993. [Google Scholar]
Stuth, J.W.; Sheffield, W.J. Determining carrying capacity for combinations of livestock, white-tailed deer and exotic ungulates. In Wildlife Managagement Handbook; Texas A&M University: College Station, TX, USA, 2001; pp. 5–12. [Google Scholar]
Stuth, J.W.; Kamau, P.N. Influence of woody plant cover on dietary selection by goats in an Acacia senegal savanna of East Africa. Small Rumin. Res. 1990, 3, 211–225. [Google Scholar] [CrossRef]
Degen, A.A. Sheep and goat milk in pastoral societies. Small Rumin. Res. 2007, 68, 7–19. [Google Scholar] [CrossRef]
Haenlein, G.F. W. The Nutritional Value of Sheep Milk. Available online: http://www.smallstock.info/issues/sheepmilk.htm (accessed on 21 October 2010).
Maltz, E.; Shkolnik, A. Milk production in the desert: Lactation and water economy in the black Bedouin goat. Physiol. Zool. 1980, 53, 12–18. [Google Scholar]
Meged’, S.S.; Torkaev, A.N.; Egorov, S.V.; Storozhuk, S.I. Milk productivity of breeding sheep of the Altai finewool breed. Russ. Agric. Sci. 2008, 34, 52–54. [Google Scholar] [CrossRef]
Thomson, E.F.; Bahhady, F.; Termanini, A.; Mokbel, M. Availability of home-produced wheat, milk products and meat to sheep-owning families at the cultivated margin of the NW Syrian steppe. Ecol. Food Nutr. 1986, 19, 113–121. [Google Scholar] [CrossRef]
Harris, F. Management of manure in farming systems in semi-arid West Africa. Exp. Agric. 2002, 38, 131–148. [Google Scholar] [CrossRef]
Bernabeu Auban, J.; Molina Balaguer, L.; Orozco Köhler, T.; Díaz Castillo, A.; Barton, C.M. Early Neolithic at the Serpis Valley, Alicante, Spain. In The Early Neolithic in the Iberian Peninsula. Regional and Transregional Components, Proceedings of the XV World Congress, Lisbon, Portugal, 4–9 September 2006; Diniz, M., Ed.; BAR International Series: Oxford, UK, 2008; pp. 53–59. [Google Scholar]
Bryson, R.A.; Bryson, R.U. High resolution simulations of regional Holocene climate: North Africa and the Near East. In Third Millennium BC Climate Change and Old World Collapse; Nato ASI Series I: Global Environmental Change; Nüzhet Dalfes, H., Kukla, G., Eds.; Springer Verlag: Berlin, Germany, 1996; Volume 49, pp. 565–593. [Google Scholar]
Ruter, A.; Arzt, J.; Vavrus, S.; Bryson, R.A.; Kutzbach, J.E. Climate and environment of the subtropical and tropical Americas (NH) in the mid-Holocene: Comparison of observations with climate model simulations. Quat. Sci. Rev. 2004, 23, 663–679. [Google Scholar] [CrossRef]
Arıkan, B. Macrophysical climate modeling, economy, and social organization in Early Bronze Age Anatolia. J. Archaeol. Sci. 2014, 43, 38–54. [Google Scholar] [CrossRef]
Ullah, I. The Consequences of Human Land-Use Strategies during the PPNB-LN Transition: A Simulation Modeling Approach. Ph.D. Thesis, Arizona State University, Tempe, AZ, USA, 2013. [Google Scholar]
Bogaard, A. “Garden agriculture” and the nature of early farming in Europe and the Near East. World Archaeol. 2005, 37, 177–196. [Google Scholar] [CrossRef]
Conolly, J.; Colledge, S.; Shennan, S. Founder effect, drift, and adaptive change in domestic crop use in early Neolithic Europe. J. Archaeol. Sci. 2008, 35, 2797–2804. [Google Scholar] [CrossRef]
Colledge, S.; Conolly, J.; Shennan, S. The evolution of Neolithic farming from SW Asian origins to NW European limits. Eur. J. Archaeol. 2005, 8, 137–156. [Google Scholar] [CrossRef]
Halstead, P. Traditional and ancient rural economy in Mediterranean Europe: Plus ça change? J. Hell. Stud. 1987, 107, 77–87. [Google Scholar] [CrossRef]
Butzer, K.W. Archaeology as Human Ecology; Cambridge University Press: Cambridge, UK, 1982. [Google Scholar]
McClure, S.B.; Barton, C.M.; Jochim, M.A. Human behavioral ecology and climate change during the transition to agriculture in Valencia, eastern Spain. J. Anthropol. Res. 2009, 65, 253–269. [Google Scholar] [CrossRef]
McClure, S.; Jochim, M.A.; Barton, C.M. Behavioral ecology, domestic animals, and land use during the transition to agriculture in Valencia, eastern Spain. In Foraging Theory and the Transition to Agriculture; Kennett, D., Winterhalder, B., Eds.; Smithsonian Institution Press: Washington, DC, USA, 2006; pp. 197–216. [Google Scholar]

© 2015 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 license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Barton, C.M.; Ullah, I.; Heimsath, A. How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes. Land 2015, 4, 578-606. https://doi.org/10.3390/land4030578

AMA Style

Barton CM, Ullah I, Heimsath A. How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes. Land. 2015; 4(3):578-606. https://doi.org/10.3390/land4030578

Chicago/Turabian Style

Barton, C. Michael, Isaac Ullah, and Arjun Heimsath. 2015. "How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes" Land 4, no. 3: 578-606. https://doi.org/10.3390/land4030578

APA Style

Barton, C. M., Ullah, I., & Heimsath, A. (2015). How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes. Land, 4(3), 578-606. https://doi.org/10.3390/land4030578

Article Menu

How to Make a Barranco: Modeling Erosion and Land-Use in Mediterranean Landscapes

Abstract

1. Introduction

2. MedLanD Modeling Laboratory (MML)

2.1. Surface Process Model

2.2. Human Land-Use Model

3. Modeling Experiment Design

3.1. Creating the Digital Landscape

3.2. Modeling Past Climate

3.3. Modeling Land-Use

3.4. Sensitivity Testing and Experiment Repetition

4. Simulation Results

4.1. Erosion and Deposition in Barrancos

4.2. Effects of Land-Use on Barranco Incision

5. Conclusions

Acknowledgments

Author Contributions

Appendix

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI