**Contents**


## **About the Editor**

**Manoj K. Jha** is an associate professor in the Department of Civil, Architectural and Environmental Engineering at North Carolina A&T State University. Dr. Jha has about 20 years of experience working in multi-scale modeling of hydrology and water resources/quality, evaluating water management strategies and impact assessment studies due to urban and agriculture best-management practices, land-use change and climate variability. He has made significant contributions in several areas of water resource research nationally and internationally. Specifically, he has developed and applied various fields to watershed scale models to extend the boundaries of our knowledge and understand the impacts of land use and climate change on hydrology, water availability, and water quality. The methods and tools developed by Dr. Jha have made an impact on how we should manage our watersheds to help solve the current water quality problems. Dr. Jha has made significant scientific contributions, as shown by the number of publications in reputable journals and the numerous appearances at professional meetings and the very high and respectable citation index. He has served as a PI/Co-PI on several research grants and contracts worth millions of dollars in funding. His high level of productivity in scientific contributions is reflected by his many awards and honors, including the recent 2019 Senior Researcher of the Year Award in 2019.

## **Preface to "Impacts of Landscape Change on Water Resources"**

In order to manage our valuable water resources, it is imperative to understand the degree of vulnerabilities and resiliency towards changes in the landscape. Continuous changes in land use and land cover can have many drivers, including population growth, urbanization, demand for food, evolution of socio-economic structure, policy regulations, and climate variability. Potential impacts due to these changes could range from changes in water availability (due to changes in losses of water to evapotranspiration and recharge) to degradation in water quality (increased erosion, salinity, chemical loadings, and pathogens). Fields studies are conducted to understand this complexity at local scales, while analyses at regional or watershed scales adopt modeling and simulation strategies. A range of tools, including hydrological, biophysical and ecosystem models, are used (stand-alone or in combination) to investigate important questions regarding impacts in order to inform the decision-making process. These decision analysis tools identify landscape-change impacts, risks, and uncertainties to provide guidance in making key management decisions. In this Special Issue, we include research and discussion topics from field investigations, as well as analytical and modeling studies to better understand the connection between landscape change and water resources at various scales.

> **Manoj K. Jha** *Editor*

### *Editorial* **Impacts of Landscape Changes on Water Resources**

#### **Manoj K. Jha**

Civil, Architectural, and Environmental Engineering, North Carolina A&T State University, Greensboro, NC 27411, USA; mkjha@ncat.edu

Received: 30 July 2020; Accepted: 7 August 2020; Published: 10 August 2020

**Abstract:** Changes in land use and land cover can have many drivers, including population growth, urbanization, agriculture, demand for food, evolution of socio-economic structure, policy regulations, and climate variability. The impacts of these changes on water resources range from changes in water availability (due to changes in losses of water to evapotranspiration and recharge) to degradation of water quality (increased erosion, salinity, chemical loadings, and pathogens). The impacts are manifested through complex hydro-bio-geo-climate characteristics, which underscore the need for integrated scientific approaches to understand the impacts of landscape change on water resources. Several techniques, such as field studies, long-term monitoring, remote sensing technologies, and advanced modeling studies have been contributing to better understanding the modes and mechanisms by which landscape changes impact water resources. Such research studies can help unlock the complex interconnected influences of landscape on water resources for quantity and quality at multiple spatial and temporal scales. In this Special Issue, we published a set of eight peer-reviewed articles elaborating on some of the specific topics of landscape changes and associated impacts on water resources.

**Keywords:** landscape change; water resources analysis; water modeling; impact assessment

#### **1. Introduction**

Landscape change and its impact on water resources is a vast topic that encompasses fundamental and applied research in multiple dimensions including water resources science and engineering, agriculture, geology, geography, economics, and social sciences. Landscape change can have many manifestations, such as changes due to urbanization, industrialization, commercialization of marginal lands, agriculture, farmers' decisions on the use of croplands, increased use of land through deforestation and drainage of wetlands, government policy decisions for environmental regulations, natural disasters such as floods and droughts, changing climatic and environmental conditions, and others. Subsequently, the impacts of these changes in one form or another on water resources are realized at spatial scales (local impacts contributing to regional scales), temporal scales (short-term vs. long-term changes), changes in water footprints through changes in hydrological processes, and changes in water/environmental quality (sediments, nutrients, and pathogens). For example, changing land use through deforestation and/or drained wetland will influence changes in infiltration and runoff characteristics, thereby affecting evapotranspiration, groundwater recharge, and sediment and water yield [1–4].

It is well-known that the changes in land use have large impacts on water resources; however, quantifying these impacts remains among the more challenging problems in managing water resources [5]. One of the major challenges is the complex interconnection of water within the hydro-bio-geo-climate characteristics [6]. As land use changes, it will alter the water balance through changes in groundwater recharge, runoff, and evapotranspiration. Water movement will be affected due to changes in soil physical properties such as moisture content and soil temperature. Variety of land use types will have associated changes in land characteristics affecting the movement of water in

variety of ways. Added to this challenge is the response time of the impacts. For example, groundwater systems response to changes in land use may vary widely from days to decades. Similarly, large scale changes in land use will impact evapotranspiration to a large extent that will be propagated through hydrologic systems and may potentially modify regional weather patterns and climate variability in an unknown future time.

Technological advancements over the years and continuous research efforts have pushed the boundary of science to better understand and quantitatively assess the impacts of landscape change on water resources. While field studies have been proven useful to understand the complexity of impacts at local scales [7], analyses at regional or watershed scales adopt modeling and simulation strategies [8,9]. A range of tools, including hydrological, biophysical, ecosystem models have been developed and used (stand alone or in combination) for investigation and inform the decision-making process. These decision analysis tools identify landscape-change impacts, risks, and uncertainties to provide guidance to make key management decisions [10,11].

This Special Issue presents studies [12–19] that describe the application of a variety of observational and modeling tools and techniques to evaluate the impacts of landscape changes on water resources in watersheds. Landscape changes included change due to agricultural best management practices (BMPs), low impact development (LID) in urban settings, conservation agriculture practices, conversion of erosion hot spot cropland into forest, and others. The impact on water resources included the changes in streamflow, stream and soil temperature, evapotranspiration, sediment, and others. The application of methods included study watersheds in U.S., Hungary, China, South Korea, and Ethiopia.

#### **2. Summary of Papers in the Special Issue**

Table 1 shows the comparative analysis of all articles of this Special Issue in terms of the type of the land use change analysis, evaluated impact assessment parameters, and the methods used. Brief summaries of each of the eight published papers are also presented.




**Table 1.** *Cont*.

#### *2.1. Cumulative E*ff*ects of Low Impact Development on Watershed Hydrology in a Mixed Land-Cover System, by Hoghooghi et al., 2018*

LID practices are designed to reduce the impact of land use change on hydrology. This study used a spatially explicit ecohydrological model VELMA, to assess the impact of LID techniques at the watershed scale. The authors calibrated and validated the model for streamflow. Hydrological effects of three common LID practices (rain gardens, permeable pavement, and riparian buffers) were tested on a 0.94 km2 mixed land cover semi-urban watershed (27 percent impervious) in Ohio, USA. LID practices were shows to perform as expected for effectively reducing the peak flow and increasing infiltration but with limited efficiency in semi-urban watershed.

#### *2.2. Modeling Landscape Change E*ff*ects on Stream Temperature Using the Soil and Water Assessment Tool, by Mustafa et al., 2018*

Stream temperature is an important factor in regulating fish behavior and habitat. This study investigates the impact of landscape change on stream temperature. The authors developed a mechanistic stream temperature module within the watershed modeling environment of the SWAT model and applied it to the 782 km<sup>2</sup> watershed in Oregon, USA. The model was calibrated for flow and stream temperature before examining the changes in stream temperature due to change in land use. The model was able to capture the increased stream temperatures in agricultural sub-basins compared with forested sub-basins.

#### *2.3. Improved Soil Temperature Modeling Using Spatially Explicit Solar Energy Drivers, by Halama et al., 2018*

Soil temperature affects ecosystem properties including increasing the water temperature. This study demonstrated that local solar energy information improved soil temperature modeling estimates simulated by a soil temperature subroutine within a larger ecohydrological watershed model VELMA. Authors calibrated the model using the data available from Oregon's Crest-to-Coast Environmental Monitoring Transect (O'CCMoN) sites. Results demonstrated the benefit of including spatially explicit representations of solar energy within watershed-scale models that simulate soil temperature.

#### *2.4. Issues of Meander Development: Land Degradation or Ecological Value? The Example of the Sajó River, Hungary, by Bertalan et al., 2018*

River channels and their surrounding floodplains enhance landscape evolution and the diversification of environments. This study investigates the geomorphological development and effects of bank erosion along meandering Sajó River in Hungary. Authors performed GIS analysis of three consecutive meandering bends over 10 periods between 1952 and 2017 based on archive aerial imagery and UAV-surveys. Analyses revealed that the meandering (channel sinuosity) was directly proportional to forest density (dominant, compact, and connected) which provided high ecological value.

#### *2.5. E*ff*ects of Di*ff*erent Spatial Configuration Units for the Spatial Optimization of Watershed Best Management Practice Scenarios, by Zhu et al., 2019*

Variation in spatial configurations of BMPs at the watershed scale may have significantly different environmental effectiveness. This study investigated and compared the effects of four main types of spatial configuration units for BMP scenarios optimization. Optimization was conducted based on a fully distributed watershed modeling framework, the Spatially Explicit Integrated Modeling System (SEIMS) using SWAT, and an intelligent optimization algorithm Non-dominated Sorting Genetic Algorithm II (NSGA-II). Results showed that the different BMP configuration yielded significant differences in near-optimal Pareto solutions, optimizing efficiency, and spatial distribution of BMP scenarios. BMP configuration units that support the adoption of expert knowledge on the spatial relationships between BMPs and spatial locations (e.g., hydrologically connected fields, slope position units) are considered to be the most valuable spatial configuration units for watershed BMP scenarios optimization and integrated watershed management.

#### *2.6. Evaluating the E*ff*ectiveness of Spatially Reconfiguring Erosion Hot Spots to Reduce Stream Sediment Load in an Upland Agricultural Catchment of South Korea, by Choi et al., 2019*

Soil erosion has a negative impact on the environment and socioeconomic factors by degrading the quality of both nutrient-rich surface soil and water. This modeling study demonstrated the effectiveness of converting soil erosion hot spots within the watershed into forest for reducing the sediment yield significantly.

#### *2.7. Scaling-Up Conservation Agriculture Production System with Drip Irrigation by Integrating MCE Technique and the APEX Model, by Assefa et al., 2019*

Conservation agriculture, which promotes no-till, mulching, and diverse cropping, provides higher water use efficiency in addition to improving soil fertility and crop yield. This study demonstrated the scaling-up impacts of conservation agriculture on a regional scale. The calibrated biophysical model APEX in combination with GIS-based multi-criteria evaluation (MCE) technique was used to extend the modeling analysis to the national scale in Ethiopia. Results indicated that the conservation agriculture with drip irrigation technology could improve groundwater potential for irrigation up to five folds and intensify crop productivity by up to three to four folds across the nation.

#### *2.8. Flooding Urban Landscapes: Analysis Using Combined Hydrodynamic and Hydrologic Modeling Approaches, by Jha and Afreen, 2020*

Urban landscape dictates the extent of inundation during a flood event, affecting vulnerable infrastructures. This study presents a systematic approach of combining hydrodynamic model HEC-RAS with hydrologic model SWAT in delineating flood inundation zones, and subsequently assessing the vulnerability of critical infrastructures in the Blue River Watershed in Kansas City, Missouri. Results demonstrate the usefulness of such combined modeling systems to predict the extent of flood inundation and thus support analyses of management strategies to deal with the risks associated with critical infrastructures in an urban setting.

#### **3. Conclusions**

Landscape changes have direct linkages with changes in hydrology in terms of water balance components. As land use characteristics change, it will alter the hydrology at the local scale, leading to the impacts on water availability and associated water quality to regional scales and at various temporal scales. The papers in this Special Issue describe the applications of a variety of observational and modeling tools and techniques to evaluate the impacts of landscape changes on water resources. The studies can be categorized into four subject areas: (1) impact assessment due to implementation of management practices [12–15], (2) impact of landscape on stream and soil temperature [16,17], (3) landscape and river meandering [18], and (4) landscape for flood inundation [19].

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

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


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

## *Article* **Cumulative Effects of Low Impact Development on Watershed Hydrology in a Mixed Land-Cover System**

**Nahal Hoghooghi 1,2, Heather E. Golden 3,\*, Brian P. Bledsoe 2, Bradley L. Barnhart 4, Allen F. Brookes 4, Kevin S. Djang 5, Jonathan J. Halama 4, Robert B. McKane 4, Christopher T. Nietch <sup>6</sup> and Paul P. Pettus <sup>4</sup>**


Received: 20 June 2018; Accepted: 20 July 2018; Published: 27 July 2018

**Abstract:** Low Impact Development (LID) is an alternative to conventional urban stormwater management practices, which aims at mitigating the impacts of urbanization on water quantity and quality. Plot and local scale studies provide evidence of LID effectiveness; however, little is known about the overall watershed scale influence of LID practices. This is particularly true in watersheds with a land cover that is more diverse than that of urban or suburban classifications alone. We address this watershed-scale gap by assessing the effects of three common LID practices (rain gardens, permeable pavement, and riparian buffers) on the hydrology of a 0.94 km<sup>2</sup> mixed land cover watershed. We used a spatially-explicit ecohydrological model, called Visualizing Ecosystems for Land Management Assessments (VELMA), to compare changes in watershed hydrologic responses before and after the implementation of LID practices. For the LID scenarios, we examined different spatial configurations, using 25%, 50%, 75% and 100% implementation extents, to convert sidewalks into rain gardens, and parking lots and driveways into permeable pavement. We further applied 20 m and 40 m riparian buffers along streams that were adjacent to agricultural land cover. The results showed overall increases in shallow subsurface runoff and infiltration, as well as evapotranspiration, and decreases in peak flows and surface runoff across all types and configurations of LID. Among individual LID practices, rain gardens had the greatest influence on each component of the overall watershed water balance. As anticipated, the combination of LID practices at the highest implementation level resulted in the most substantial changes to the overall watershed hydrology. It is notable that all hydrological changes from the LID implementation, ranging from 0.01 to 0.06 km<sup>2</sup> across the study watershed, were modest, which suggests a potentially limited efficacy of LID practices in mixed land cover watersheds.

**Keywords:** LID practices; watershed scale; impervious area; peak flow; surface runoff; shallow subsurface runoff and infiltration; evapotranspiration

#### **1. Introduction**

Urbanization alters natural hydrological systems by altering stream channel networks (e.g., channelization and burial), creating microclimates (e.g., urban heat islands), and generating rapid runoff from precipitation and snowmelt events [1]. These changes have direct impacts on surface and groundwater quantity and quality. Conventional urban stormwater management practices are often developed to control runoff and minimize flooding; however, these systems can be costly and may not directly address issues, such as reductions in infiltration and groundwater storage via impervious surfaces that may lead to urban flooding, erosion, and the degradation of water quality [2].

In recent years, alternative stormwater management practices, such as Low Impact Development (LID), have been adopted (e.g., bioretention cells or rain gardens, permeable pavements, and bioswales) in many urban and suburban areas [3,4]. LID, also called sustainable urban drainage systems (SUDSs), among other globally varying names [5], is an approach that uses soils, vegetation, and landscape design to control nonpoint source runoff and pollutants in urban systems. A goal of LID is to promote watershed resilience through "green" design [6].

There is a growing body of literature focused on evaluating the local (e.g., plot or site) scale effectiveness of LID. Several recent papers have synthesized the key findings of studies assessing the effects of different LID practices, including field experiments and modeling studies, on water quantity (e.g., peak flow and runoff volume) and water quality (e.g., nitrogen, phosphorous, and total suspended solids) at local scales [7–11]. These previous studies provide foundational research for scaling LID approaches to watersheds. However, limited evidence of LID effectiveness at the watershed scale exists [12,13], and research focusing on LID impacts at watershed scales is just beginning to emerge [14]. Therefore, questions remain about how LID practices can individually or cumulatively affect watershed hydrology [14,15].

Experimental studies, designed to investigate the watershed-scale effects of LID, have provided critical insights into how watershed hydrology responds to these approaches. For example, Jarden et al. [16] designed a paired watershed approach to quantify the effect of street-connected bioretention cells, rain gardens, and rain barrels on peak discharge and total storm runoff. The results from the subwatershed with smaller LID lots and underdrain connections showed a substantial reduction in peak discharge (up to 33%) and total storm runoff (up to 40%). Additionally, recent field-based research provides evidence of the cumulative watershed scale effects of LID on hydrologic responses, such as peak flows and pollutant loads [17–20]. Such experimental studies can be resource intensive (e.g., financial, personnel, time) [21]; however, process-based models provide a means to go beyond measured data and explore the projected "what if" LID scenarios using potentially less resources.

Process-based or mechanistic watershed models, which simulate hydrological (and other) processes and outputs for different water balance components (e.g., streamflow, evapotranspiration), are critical tools to understand the influence of LID practices on watershed processes [22]. The Storm Water Management Model (SWMM) [23] has been used to simulate the effects of different LID practice implementations (porous pavement, rain barrels, and rain garden) on runoff and flood risk reductions in an 87.6 km<sup>2</sup> urban watershed and model the cumulative effects of street-side bioretention cells, rain gardens, and rain barrels in a 0.12 km2 residential watershed [21]. Overall, the results indicate increases in evaporation and infiltration, as well as decreases in surface runoff and discharge, across different return periods. The performance of LID practices (rain gardens, permeable pavements, and rainwater harvesting tanks) has also been evaluated under different urban land use densities using the Soil and Water Assessment Tool (SWAT) [24], demonstrating that the effectiveness of LID practices differs among the urban land use densities [25].

The aforementioned studies and others (e.g., [26–29]) advance current knowledge on the effectiveness of LID practices at watershed scales using various process-based model approaches; however, all of these studies focus on watersheds that are entirely urban or suburban. A clear need exists for an understanding of the extent to which LID approaches are effective in mixed land cover watersheds, i.e., those with urban and suburban land cover *in addition to others* (e.g., forest

and agriculture). Furthermore, a spatially explicit approach toward representing LID practices and associated hydrological processes to analyze the effects of varying patterns of mixed land use and land cover under different management practices is critical [30], as most approaches are challenged with representing spatial landscape heterogeneities [31].

In this paper, we assess how LID implementation affects watershed hydrologic responses in a mixed land cover watershed. Specifically, we ask: How does the type and extent of LID practices affect water balance components, including surface runoff, peak flows, evapotranspiration, shallow subsurface flow, and infiltration, in a mixed land cover watershed? We do this by using a spatiallyexplicit ecohydrological model, called Visualizing Ecosystems for Land Management Assessments (VELMA) [32] for a variety of scenarios associated with LID and the implementation of forested riparian buffers. Our study is one of the first, to our knowledge, to examine LID implementation at the watershed scale using spatially explicit modeling approaches in a system with mixed suburban, agricultural and forest land cover. As a result, we discuss the implications of this study for effective stormwater management in mixed land cover systems and future research directions toward this goal.

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

#### *2.1. Study Area Description*

The Shayler Crossing (SHC) watershed is a subwatershed of the East Fork Little Miami River Watershed in southwest Ohio, USA and falls within the Till Plains region of the Central Lowland physiographic province. The Till Plains region is a topographically young and extensive flat plain, with many areas remaining undissected by even the smallest stream. The bedrock is buried under a mantle of glacial drift 3–15 m thick [33,34]. The Digital Elevation Model (DEM) has a maximum value of ~269 m (North American\_1983 datum) within the watershed boundary (Figure 1). The soils are primarily the Avonburg and Rossmoyne series, with high silty clay loam content and poor to moderate infiltration [35]. Average annual precipitation for the period, 1990 through 2011, was 1097.4 ± 173.5 mm. Average annual air temperature for the same period was 12 ◦C [36].

We considered SHC a mixed land cover watershed, located on the east side of Cincinnati, Ohio, with a drainage area of 0.92 km<sup>2</sup> (Figure 1). The primary land uses consist of 64.1% urban or developed area (including 37% lawn, 12% building, 6.5% street, 6.4% sidewalk, and 2.1% parking lot and driveway), 23% agriculture, and 13% deciduous forest (Table 1). Total imperviousness covers approximately 27% of the watershed area, the majority of which is directly connected to a storm sewer system without any intermediary controls [30]. The watershed was chosen for this study because it is part of the East Fork Little Miami River Watershed, where a long-term monitoring and focused modeling effort is being conducted by the US Environmental Protection Agency (EPA), Office of Research and Development (ORD), Ohio Environmental Protection Agency (Ohio EPA), and Clermont County (Ohio) Stormwater Division.

#### *2.2. Input Data*

We obtained average daily precipitation and temperature data from a weather station, located approximately 13 km from the north boundary of the watershed at 84.2909◦ W, 39.194◦ N [37]. Streamflow has been monitored, from 3 April 2006 to the present day, using a stage sensor (600 LS Sonde with temperature, conductivity, and shallow vented level sensors, YSI Inc., Yellow Springs, OH, USA) at the watershed outlet. Water depth was recorded at 10-min intervals and converted to streamflow (m3 s<sup>−</sup>1) using a rating curve, developed by US EPA. We obtained a 10 m resolution DEM, Soil Survey Spatial Tabular (SSURGO 2.2) soil data, and National Land Cover Dataset (NLCD) land use data from the Natural Resources Conservation Service (NRCS) Geospatial Data Gateway [38]. We further used an impervious area shape file from Clermont County, Ohio through the Center for Urban Green Infrastructure Engineering, Inc. (Milford, OH, USA).

**Figure 1.** The study area is a 0.92 km<sup>2</sup> subwatershed (Shayler Crossing) of the East Fork Little Miami River watershed, located on the east side of Cincinnati, Ohio in Clermont County, USA. The watershed outlet is identified as a red triangle.


**Table 1.** Summary of Shayler Crossing watershed land use area and characteristics.

#### *2.3. Model Description*

To simulate the effect of LID on watershed hydrology, we used the Visualizing Ecosystems for Land Management Assessments (VELMA) model. VELMA is a spatially distributed ecohydrological model that couples watershed hydrology and carbon (C) and nitrogen (N) cycling in plants and soils, and the transport of water, C, and N from the terrestrial landscape to streams [32]. VELMA is not an "urban hydrology" model according to the strict tradition of stormwater management models (e.g., SWMM). Its key strengths are its spatially explicit representation of hydrological and biogeochemical processes and broad applicability to a variety of ecosystems, such as forest, agricultural, and urban, in order to assess the effects of LID in mixed land cover systems. Urban LID practices can be represented in the model using modifications to present watershed permeability, lateral and horizontal hydraulic conductivities, and land cover (see Section 2.5). VELMA's spatially explicit

grid-based structure affords the capacity to represent transitions from directly connected to indirectly connected impervious areas by replacing values on a cell by the cell basis for the aforementioned model representations. The model is also capable of scaling hydrologic and biogeochemistry responses across multiple spatial (hillslopes to basins) and temporal (days to centuries) scales [21]. VELMA's visualization and interactivity features are packaged in an open-source, open-platform programming environment (Java/Eclipse) [32].

VELMA's modeling domain is a three-dimensional matrix that includes information regarding surface topography, land use, and four soil layers. VELMA uses a distributed soil column framework to model the lateral and vertical movement of water and nutrients through the four soil layers. A soil water balance is solved for each layer. The soil column model has three coupled submodels: (1) A hydrological model that simulates the vertical and lateral movement of water within the soil and losses of water from soil and vegetation in the atmosphere; (2) a soil temperature model that simulates daily soil layer temperatures based on surface air temperature; and (3) a biogeochemistry model that simulates C and N dynamics.

A simple logistical function, based on the degree of saturation, is applied to capture the breakthrough characteristic of soil water. Potential evapotranspiration (PET) is estimated using the simple temperature-based method of Hamon [39]. Evapotranspiration (ET) increases exponentially as soil water storage increases, and it reaches the PET rate as the soil water storage reaches saturation. The VELMA simulator engine allows for the specification of a spatial data map, with permeability fractions for each grid cell value (here, each 10 m grid cell). The grid's permeability fractions are taken into account when determining how much of a cell's total water inflow (e.g., from rain, snow melt, and lateral surface movement) penetrates into the first layer of the soil column. A permeability of 0 is completely impermeable (no water penetrates from the surface to the first soil layer), and 1 is completely permeable (all water penetrates from the surface to the first soil layer).

The soil column model is placed within a watershed framework to create a spatially distributed model applicable to watersheds (Figure 2, shown here with LID practices). Adjacent soil columns interact through down-gradient water transport. Water entering each pixel (via precipitation or flow from an adjacent pixel) can either first infiltrate into the implemented LID and the top soil layer, and then to the downslope pixel, or continue its downslope movement as the lateral surface flow. Surface and subsurface lateral flow are routed using a multiple flow direction method, as described in Abdelnour et al. [21]. A detailed description of the processes and equations can be found in McKane et al. [32], Abdelnour et al. [21], Abdelnour et al. [40].

#### *2.4. Watershed Model Setup*

We used VELMA's pre-processor tool (called Java Processing Digital Elevation Model (JPDEM) [32,41]) to fill sinks, determine flow direction, and compute the flow contribution area of a 10 m DEM [32]. The watershed boundary was delineated, and the watershed outlet was assigned using VELMA's pre-processor [32]. All DEM, soil, and land use maps were clipped so that they have the same number of columns and rows for the American Standard Code for Information Interchange (ASCII grid, Esri, Inc., Esri grid format ArcGIS Desktop 10.0 Help, http://desktop.arcgis.com/en/ arcmap/10.3/manage-data/raster-and-images/esri-grid-format.htm) input in VELMA. The soil and the land use maps contained ID numbers for every cell in the simulation area, which corresponded to one or more of VELMA's simulator configurations. We assigned two of VELMA's soil configurations to represent Rossymoyne and Avonburg soil types and seven land use configurations to represent agriculture, forest, lawn, buildings, streets, sidewalks, parking lots and driveways. We merged wet pond pixels with lawn pixels because currently lakes and ponds are not implementable in VELMA.

**Figure 2.** Generalized structure of the VELMA model, and applications for LID. VELMA's domain is a watershed (**a**). Each grid cell within the watershed has four soil layers (**b**), including LID applied to the land surface). Infiltration from LID is transported to the first soil layer (black arrow). Further vertical transport of water (in this paper), carbon, and nitrogen transport can occur between each grid cell's four layers (white arrows), and surplus water (or carbon, nitrogen) is transferred from a grid cell to the adjacent, most down gradient cell(s) in the watershed (blue arrows). P = Precipitation; ET = evapotranspiration. Modified from Abdelnour et al. [21].

#### *2.5. Base Model Parameterization, Calibration and Validation*

We performed the base model calibration with daily streamflow at the outlet of the watershed from 1 January 2009 to 31 December 2010 with 2008 as a model warmup period and from 1 January 2011 to 31 December 2011 as a validation period. Our calibration period (2009 and 2010) included normal precipitation years (1040 and 1046 mm), and our verification period (2011) was a wet year (1660 mm). We defined a 'wet' period as greater than one standard deviation from the mean precipitation (>1270.1 mm) and a dry period as less than one standard deviation (<923.9 mm).

Calibration was conducted through both semi-automatic and manual calibrations. We used autocalibration to screen for sensitive parameters and reduce the solution space. Manual calibration was implemented as a second phase to further refine the parameter values. For the initial automatic calibration, we used the MOEA-VELMA calibration tool that links VELMA with the Multiobjective Evolutionary Algorithm (MOEA) [42] framework in Java. The MOEA framework uses evolutionary algorithms to solve multiobjective optimization problems, and the MOEA-VELMA calibration tool leverages this ability to tune model input parameters to minimize the differences between simulated results and observed data. Several parameters were chosen to calibrate the model, including soil layer thickness, saturated hydraulic conductivity, porosity fraction, bulk density, wilting point, field capacity, and PET parameters. The MOEA-VELMA calibration tool then implemented NSGA-II [43], using the MOEA framework, and searched for the optimal set of input parameters to optimize our objective function, that is, Nash Sutcliffe Efficiency (NSE) [44] for the observed and predicted daily streamflow:

$$NSE = 1 - \frac{\sum\_{i=1}^{n} \left| O\_i - P\_i \right|^2}{\sum\_{i=1}^{n} \left| O\_i - \overline{O} \right|^2} \tag{1}$$

where *Oi* is the ith measured variable (e.g., discharge), *Pi* is the ith predicted variable, *O* is the arithmetic average of the measured variable, and *n* is the total number of observations. The NSE coefficient ranges between 1 (perfect fit) and negative infinity. An efficiency below zero implies that the mean value of the observed value is a better predictor than the model.

After almost 500 simulations, we narrowed the range of selected sensitive parameters and ran the MOEA-VELMA calibration tool for an additional 500 simulations. Then, we picked the solutions with a higher NSE and used those parameter ranges in the manual calibration.

After the initial semi-automatic calibration, we conducted manual calibration through visual analysis to capture trends in observed streamflow, using NSE in addition to percent bias (*PBIAS*) [45] and root mean squared error (RMSE) [46]. PBIAS measures the average tendency of the predicted data to be larger or smaller than observed values. It is also measures over- and underestimation of bias [44]:

$$PBIAS = \frac{\sum\_{i=1}^{n} (O\_i - P\_i)}{\sum\_{i=1}^{n} O\_i} \times 100\tag{2}$$

and RMSE is the square root of the mean square error and varies from zero to large positive values:

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(O\_i - P\_i\right)^2} \tag{3}$$

To ensure the simulations provided reasonable volumetric matches with observed data, we also used a total simulated to total observed annual streamflow ratio (Sim:Obs) for each simulation year. If the Sim:Obs was >1, simulated streamflow from the year exceeded that of the observed streamflow. If it was <1, the opposite was true, and Sim:Obs = 1 suggested a perfect match between the total annual simulated and observed streamflow.

The soil thickness of each layer was parameterized using United States Department of Agriculture (USDA) soil survey data for the study area [35]. We used the MOEA-VELMA calibration tool to calibrate saturated vertical and horizontal hydraulic conductivities for each soil layer. Other soil physical characteristics (porosity, field capacity, wilting point, and bulk density) were obtained based on soil texture class (Table 2) [32]. We obtained the first term of the PET *Hamon* equation (petParam1) for different cover types using the MOEA-VELMA calibration tool, with the second term of *Hamon* equation set to a constant value of 0.622, based on Abdelnour et al. [21]. A *be* parameter is a calibration constant; it is an ET coefficient used in the logistic equation that computes ET from PET. We estimated this parameter value from autocalibration. Air density (*roair*) was constant and set to 1300 g m−<sup>3</sup> (Table 3). We adjusted all soil physical characteristics and PET parameters to best match the observed streamflow during manual calibration. The parameters and their final model values are shown in Tables 2 and 3. Setting soil parameters to zero produces an error in the VELMA output; therefore, we set soil parameter values for impervious areas to those of the clay soil texture class, using the approach by McKane et al. [32].


**Table 2.** Soil parameters for the base model (\* = calibrated; all other values from McKane et al. [32]). *R* stands for Rossmoyne and *A* stands for Avonburg soil type.


**Table 3.** Calibrated Potential Evapotranspiration parameters for the base model.

#### *2.6. Low Impact Development (LID) Configurations, Scenarios, and Model Parameters*

To evaluate the effectiveness of LID practices based on the relative daily changes in watershed hydrology compared to the calibrated base model, we simulated three types of LID scenarios: Rain gardens (RG), permeable pavements (PP), and forested riparian buffers (RB). Our goal was to derive a relative understanding of how different spatial distributions of select LID types may affect hydrology in this mixed land over system; therefore, we did not aim to represent specific stakeholder-selected LID practices for the watershed (e.g., exact sites where landowners would agree to implementation).

To implement the LID scenarios, we replaced grid cells in the calibrated base model, identified as impervious, areas with one of two LID practices: RG or PP, depending on the impervious area type (see below). We further replaced grid cells in agricultural land cover along a stream with RB. We ran each scenario as a separate model using evenly distributed spatial configurations of 25%, 50%, 75% and 100% conversions for RG (in sidewalk locations) and PP (at parking lots and driveways; see Figure 3 for an example). Each spatial distribution of RG and PP met or exceeded the watershed's water quality volume for bioretention (i.e., generally speaking, the volume of water treated by LID practices to control in low to medium magnitude storm events), as recommended by the Ohio Department of Natural Resources [47]. We also placed RP at 20 m and 40 m on each side of the stream in the agricultural land of the Northern part of the watershed (Figure 1). This resulted in 10 simulated LID scenarios (4 RG, 4 PP, and 2 RB) for comparison. We note here that a large-scale conversion of impervious areas to LIDs (e.g., our 100% conversion scenarios) may not be reasonable, in terms of both financial cost and the willingness of the community [46]; however, these conversion configurations can provide a maximum mitigation potential for decision support.

RG and PP were chosen for the scenarios because they are reasonable retrofitting measures for the studied watershed, are the most promising LID practices for reductions in peak flow and runoff volume [8,48], and can be applied and assessed in the VELMA model. RBs were selected because they currently do not exist in the agricultural land of the watershed (and therefore the base model). Their addition was used for comparisons of the watershed-scale hydrological responses of RG and RB conversions on impermeable areas. We selected 20 m and 40 m buffers to go beyond Ohio EPA's

requirement that forested area must be maintained for a minimum of the first 15 m of the area on either bank [49].

We assessed the effect of the LID scenario implementation: (a) Individually and (b) using an LID combination scenario (i.e., fully implementing RG, PP, and RB with the maximum level of implementation). The individual model scenario runs of land cover conversions to RG and PP, for each spatial configuration, included: (a) Sidewalks were converted to RG and (b) parking lots and driveways were converted to PP (Table 4 and Figure 3). Lawns were not converted to RG. The percentage of the watershed that was converted to RG, PP, and RB practices at different implementation levels is shown in Table 5.


**Table 4.** LID configurations in the model and conversion levels.

**Figure 3.** The spatial configuration of different land covers in the SHC watershed, with 10m cell resolution: (**a**) Current land cover, and (**b**) after LID implementation: Conversion of a 100% spatial configuration of the sidewalks to rain gardens, and driveways and parking lots into permeable pavements, and the implementation of 40 m forest buffers along both sides of the stream on the agricultural land. The watershed outlet is identified as a red triangle.


**Table 5.** The percent of the watershed converted to each LID practice and implementation level.

To implement LID into each scenario, we parameterized the soil texture, soil physical characteristics, and PET parameter values for LID practices, based on Ohio EPA requirements (Table 6) [49]. To do this for the RG scenarios, we created soil maps with a new soil class, "RG," for each spatial configuration (i.e., 25% to 100%). The RG soil maps, one for each implementation level, replaced the sidewalk pixels of the original soil map. Soils in the new RG maps were adjusted for soil depths, texture classes, and physical parameters to represent soils associated with rain gardens. The RG soil maps were based on Ohio EPA requirements for rain gardens (Table 6), which suggest that the soil media depths of a rain garden are 60–100 cm deep with loamy sand [49]. In the updated model configurations for each implementation scenario, we assumed no underdrain pipes and no outlet pipes, which are currently not implementable in VELMA.

In addition to the soil maps, we created new land cover maps for each spatial implementation level of RG (25% to 100%). We defined a new land cover, "RG," where existing sidewalks were located. For example, at the 50% RG implementation level, 50% of sidewalk's pixels of original were defined as "RG" land cover. For each new "RG" map, we parameterized the PET parameters of "RG" land cover to lawn values (Table 7).

For PP scenarios, we generated soil maps with a new soil class, "PP" which replaced parking lots and driveways at each conversion level (25% to 100%). We modified the original soil depths, soil texture classes, and soil physical parameter values for the "PP" soil class (Table 6) using the same values at each conversation level. According to the Ohio EPA Stormwater Management Practices manual, the recommended thickness of a PP system is 40–76 cm, depending on frost depth [49]. Therefore, for PP, we parameterized the hydraulic conductivity and other soil physical parameter values of the first 100 cm of the "PP" soil class [32]. We assumed that permeable pavement is a continuous pavement system (gravel) and well maintained with no clogging issues.

For RB, we created soil maps with a new soil class (Table 6) and a new land cover class (Table 7), "RB" to replace current soils and land cover at 20 m and 40 m on each side of streams where agriculture exists. Because most riparian buffers for Ohio streams are forested [49], we parameterized the soil parameters and PET values of the buffer area in the new soil and cover maps to reflect the effect of a forest rooting system and forest canopy on infiltration and ET (Tables 6 and 7).

For RG, PP and RB, new permeability fraction maps were also created for each implementation level to replace the original permeability fraction map in each model scenario. In the new permeability fraction maps, permeability fractions of 0 for impervious surfaces, such as sidewalks, parking lots and driveways, were changed to 1 and 0.95 for RG and PP, respectively. The permeability fraction for the RB scenario was changed from 0.95 for agriculture to 1 for RB forested land cover.

Once the base model was calibrated, we ran the model for each of the 10 scenarios under the different LID spatial configurations to evaluate changes in peak flows, surface runoff, ET, subsurface runoff and infiltration, and compared them to that of the base model (existing conditions).


**Table 6.** Soil parameter values used in the LID practice scenarios. RG: Rain Garden, PP: Permeable pavement, and RB: Riparian Buffer.

<sup>1</sup> Ohio EPA [49]; <sup>2</sup> McKane et al. [32], based on the Ohio EPA recommended soil texture class of loamy sand and sand for RG and PP [49]. For the RB scenario, the values were set to the loamy sand soil texture class to represent a forest rooting system [32].



<sup>1</sup> The values for RG and RB were set to the calibrated values for lawn and forest land cover (Table 3). The value for PP set to minimum value for impervious area (Table 3).

#### **3. Results**

#### *3.1. Calibration and Validation of the Base Model*

Daily streamflow calibration suggests acceptable model results across the simulation period (Figure 4a–c). The NSE, *R*2, root mean square error (RMSE), and percent bias (PBIAS) for the calibration period (2009 and 2010) were 0.50, 0.53, 3.12 and −2.40, respectively. Moriasi et al. [50] recommended that an NSE ≥ 0.50 and PBIAS ≤ ±15 can be considered satisfactory for daily streamflow simulations. RMSE varies from 0 to large positive values. The lower the RMSE, the better the model fit [46]. The optimum value for PBIAS is zero, and low magnitude values indicate better simulations. Negative values indicate model overestimation [51]. While the daily model calibration is acceptable, it tends toward underestimating peak flows (Figure 4a,b). This is confirmed by a negative PBIAS; however, the magnitude is low, which means that the bias toward peak flow underestimation is minimal. Further, the Sim:Obs were 0.77, 1.10 and 0.96 (for 2009, 2010, 2011, respectively), all of which indicated that annual volumetric streamflow estimates in the base model were satisfactory.

**Figure 4.** Plots of (**a**,**b**) calibration (NSE = 0.54, *<sup>R</sup>*<sup>2</sup> = 0.53, RMSE = 3.12, and PBIAS = <sup>−</sup>2.40) and (**c**) validation (NSE = 0.40, *R*<sup>2</sup> = 0.48, RMSE = 5.27, and PBIAS = 13.84) of the VELMA model output at the watershed outlet. The model was calibrated at a daily time step from 1 January 2009 to 31 December 2010 and validated at a daily time step from 1 January 2011 to 31 December 2011.

Simulations during the validation period captured general daily streamflow patterns; however, the model fit was less satisfactory than the calibration period (NSE of 0.40, *R*<sup>2</sup> of 0.48, RMSE of 5.27, and PBIAS of 13.84). Moreover, visual inspection of the validation plot (Figure 4c) indicated that the calibrated parameters were less successful during 2011, suggesting that calibrated model simulations may have increased limitations during wet years.

The average annual water balance components of the calibrated base model across the watershed, for the simulation period (2009–2011), are shown in Table 8. Evapotranspiration accounts for about 44% of precipitation, which approximates the lower end of the Sanford and Selnick [52] estimates of the fraction of precipitation lost to evapotranspiration in Southwest Ohio, USA.

**Table 8.** Average annual water balance components of the base model for the entire watershed, from 2009–2011. Note that precipitation is lower than the sum of the other water balance components because of the structure of VELMA's hydrological model output, which couples shallow subsurface runoff with infiltration.


#### *3.2. LID Scenarios*

We compared the simulated water balances for the three LID practices at 25%, 50%, 75% and 100%, 20 m and 40 m implementation levels and one combined LID scenario at the maximum level of implementation. Our results suggest that LID practices decreased surface runoff and peak flow, and increased ET, shallow subsurface runoff and infiltration as the LID implementation level increased (Figure 5). However, the response varied among different LID practices (Figure 5).

**Figure 5.** Percent change for watershed water balance components for three different types of LID practices at the outlet of the watershed (RG: Rain Garden, PP: Permeable Pavement, and RB: Riparian Buffer), (**a**) peak flow; (**b**) surface runoff; (**c**) ET (evapotranspiration); and (**d**) shallow subsurface runoff and infiltration. At the maximum level of implementation (100% and 40 m) RG, PP, and RB cover 6.4%, 2.1%, and 3% of the total watershed area, respectively.

Reduction in peak flows varies from about 0.5% to 5.5% among all individual LID practice scenarios, with the high reduction observed for the RG scenario at 100% and 75% implementation levels (5.5% and 4%, respectively), followed by 50% RG and 40 m RB scenarios (Figure 5a). Surface runoff decreased across all LID scenarios, with the largest reductions resulting from the RG scenario at 25%, 50%, 75% and 100% implementation levels (7%, 10.5%, 16% and 22%, respectively) (Figure 5b). PP and RB scenarios showed smaller reductions in surface runoff, ranging from 0.4 (for 25% PP) to 3.4% (for 100% PP). Reductions in surface runoff for 40 m and 20 m RB scenarios were 1.4% and 0.6%, respectively (Figure 5b). The percentage reduction in surface runoff was more than peak flows across all scenarios.

Retrofitting the baseline model with the LID increased shallow subsurface runoff and infiltration with increasing implementation levels as shown in Figure 5c,d. ET increased 2–15% for RG and RB scenarios across implementation levels, with higher increases in the 100%, 75% and 50% RG scenarios (11%, 8% and 5%, respectively). The RG scenario resulted in higher increases for both processes in comparison to other individual scenarios. Following the same trend, PP and RB scenarios increased shallow subsurface runoff and infiltration, ranging from 0.2% to 6% for different implementation levels (Figure 5d). Changes in ET for PP scenarios were negligible (Figure 5c).

Combining the three LID practices (RG, PP, and RB) at the highest implementation levels (100% for RG and PP, and 40 m for RB) resulted in the largest reductions in peak flows and surface runoff compared to individual LID implementations. The reductions in peak flow (8.5%) were modest, but considerably greater in surface runoff (26%; Figure 5a,b). The combined LID scenario resulted in the greatest increase in ET (15%), as well as a shallow subsurface runoff and infiltration (21%), in comparison with individual LID scenarios (Figure 5c,d).

The RG scenario showed the highest reduction in peak flows in comparison with PP and RB scenarios (Figure 5a). Therefore, we compared the peak flow to the percent of reduction in peak flow after RG implementation (100% scenario) during the simulation period (Figure 6). Peak flows were defined as one standard deviation above the mean simulated daily streamflow (here, 3.18 mm). We also considered the streamflow one day after we considered the peak flows to include a portion of the falling limb of the hydrograph. The percentage reduction in peak flows after RG implementation decreased exponentially with increasing peak flow conditions *R*<sup>2</sup> of 0.47; *p*-value < 0.001 (Figure 6).

**Figure 6.** Percentage reduction in peak flow after RG implementation (100% scenario) vs. peak flow (mm) during the simulation period at the watershed outlet. Peak flow was estimated as any flow that was one standard deviation above the mean simulated daily flow for the study period or flow on the day following peak flow conditions, as described (to capture part of the falling limb of the hydrograph).

#### **4. Discussion**

#### *4.1. LID Practices and Watershed-Scale Hydrological Effects*

We assess, via spatially explicit model simulations, the relative effects of different types and configurations of LID practices on watershed hydrology in a mixed land cover system. Model simulation results suggest reductions in peak flows and surface runoff, and increases in evapotranspiration and subsurface flow and infiltration, with all spatial configurations of LID at the watershed scale. This is consistent with Gagrani et al. [19], Fry and Maxwell [53], and Avellaneda et al. [54], who reported similar effects on water balance components and peak flows after the placement of different LID practices in urban watersheds, with 42–55 percent impervious surfaces and drainage areas ranging from 0.2 km2 to 12 km2.

The magnitudes of simulated water balance responses to LID placement in our watershed study were lower than other studies in strictly urban watersheds (e.g., Fry and Maxwell [53] and Avellaneda et al. [54]) and more similar to a pilot study in a small suburban watershed (1.8 km2) of Cincinnati, Ohio, where retrofitted rain gardens and rain barrels did not result in substantial runoff reductions [55]. The more limited response in our study watershed reflects, in part, the smaller extent of urban and suburban land cover compared to studies in other watersheds. Only 27 percent of our study watershed was covered in impervious surfaces, and only 31 percent of this area was converted to LID at the highest level of implementation. Therefore, our results are not completely unexpected and may point to important scale issues regarding the extent to which LID influences hydrologic regimes in mixed land use watersheds.

The simulated RB scenario did not result in a significant effect on peak flow at the watershed outlet and other water balance components. This is likely because only 3% of the total watershed area (and 13% of the watershed's agricultural land cover) was converted to RB at the highest implementation level (40 m). This indicates that the type and extent of LID practices affects cumulative watershed-scale hydrological effects [15,56].

In the mixed land cover SHC watershed, model comparisons among LID scenarios suggested that the RG was most effective across all implementation levels at reducing runoff and peak flow, and promoting ET, compared to PP and RB. RG also exerted greater control over modifying watershed water balance components, in terms of per unit area LID conversions. For example, at the 100% implementation level, RG reduced peak flows by 73% and 2%, decreased surface runoff by 50% and 86%, and increased ET by ~100%, which was 23% more than PP and RB, respectively. PP was 28% more effective at increasing shallow subsurface runoff and infiltration than RG. Recent studies point to a similar effectiveness of RG on water balance components at watershed scales [16,18]. Studies also have shown that PP can effectively mitigate surface runoff [57,58]; however, the degree of RG and PP functionality depends on the extent of the application area of LID within the watershed [59].

We found that 100% implementation of RG across the watershed was more effective at reducing peak flows during small storms than during larger ones (Figure 6). At the plot scale, Speak et al. [60] found that a green roof runoff retention significantly decreased during high rainfall events. These indicate that the effectiveness of any type of management practice, including LID, may be exponentially diminished as it loses storage capacity and becomes saturated. This counters Wadzuk et al. [61], who concluded, at the plot scale, that antecedent soil moisture conditions and "back-to-back events" are not a primary concern for biofiltration rain garden and green roof practices in recovering their infiltration capacity. However, based on our results, this finding may not be transferable to watershed scales, especially when LID practices receive both precipitation and appreciable surface runoff loading [62]. Our findings also highlight the potential importance of RG in controlling the first flush of pollutant loads and channel erosion [63,64] during more frequent storm events and a need for future research on the impact of LID practices across variable sequences of wet weather events.

#### *4.2. Implications for Stormwater Management and Future Research*

Our study provides scientists and watershed managers a glimpse into the potential influence of LID practices in mixed land cover systems, where only a portion of the watershed is converted to LID. Watershed-scale models, such as VELMA, can provide a physically-based and systematic means of projecting and evaluating the influence of various LID configurations in heterogeneous watersheds. Using this approach, our results suggest modest to minimal changes in most components of the water balance in response to LID, though these responses may be much more considerable if the watershed was exclusively urban or suburban land cover and converted to LID.

It is important to note that the location of the LID implementation with respect to the watershed outlet may also be critical [56]. For example, in our study, the watershed in agricultural areas are in the headwaters and LID implementation is in the lower portion of the watershed. Therefore, while our results suggest a limited shift in watershed hydrological dynamics with LID implementation, if LID was implemented into the upgradient of agricultural or forested land, then the magnitude of the response may be even less substantial due to attenuation from downgradient watershed processes. Based on these results, we suggest that future research needs to evaluate the hydrological effects of LID using distributed models, with a particular focus on how configurations of different land cover types influence watershed-scale LID responses, retiming of runoff delivery from subbasins of differing land cover, and antecedent soil moisture, as affected by storm sequences.

Model selection for assessing the hydrological effects of watershed-scale LID implementation is challenging because it involves trade-offs in achieving the necessary fidelity (i.e., the extent to which the model faithfully represents the modeled system) to hydrological, biological, and biogeochemical processes for prediction accuracy, while minimizing complexity and uncertainty [14]. Careful consideration of these tradeoffs is needed for future work that addresses how LID affects watershed-scale hydrological processes, particularly in mixed land cover systems. For example, existing models that explicitly integrate LID practices have been developed for urban systems and have specific LID modules for urban-based hydrological processes (e.g., SWMM and Green Infrastructure Flexible Model (GIFMod) [65]). On the other hand, models that have been explicitly developed to assess the effects of LID practices in mixed land cover watersheds (e.g., VELMA and Regional Hydro-Ecological Simulation System (RHESSys) [66]) may have a strong biogeochemical module (because of their mixed land cover focus) but a more limited hydrological capacity to physically represent LID practices, as compared with a model such as SWMM. Responses to these challenges are evolving by incorporating LID modules within ecohydrological models, such as VELMA, that provide mechanistic representations of LID performance [14] and coupling models to quantify how LID affects the fate and transport of various pollutants, as well as couple SWMM with other watershed models to improve simulations of urban hydrology [67].

Our findings suggest a clear need for the evaluation of the influence and benefits of LID in the context of other watershed land uses and their associated management [68]. For example, the incremental influence of LID on overall watershed responses, relative to management targets at different locations along the stream network, should be assessed. In this example, if a management goal is peak attenuation at the watershed outlet, a cost-benefit analysis, of how "best" to manage diffuse sources of runoff across different land cover types for peak flow reduction, would be beneficial.

Advancing the scientific understanding of the hydrological responses to LID in mixed land cover systems and linkages with the provision of diverse benefits is imperative because of the large number of watersheds globally that have mixed land cover. Future research may focus on upscaling fine scale studies to watersheds, applying a host of hydro-ecological models with LID modules or model parameter representations to address LID challenges in suburban watersheds (e.g., to provide multiple lines of evidence to support predicted outcomes), understanding LID's role in modifying baseflow (e.g., Bhaskar et al. [20]), and advancing these studies across diverse physiographic regions. Furthermore, given that we simulated and interpreted LID effects in this relatively small mixed land cover watershed, future research that applies continuous model simulations to project the hydrological

effects of different LID configurations in mixed land cover watersheds with even greater complexity than ours will help to set realistic expectations for long term LID performance in these systems. Finally, research is also needed that expands our approaches to quantifying the effects of LID practices on nutrient and sediment loads, and links LID modules within watershed-scale ecohydrology models with simulated or measured in-stream processes.

#### **5. Conclusions**

We provide one of the first studies, to our knowledge, that assesses the relative watershed-scale hydrological effects of different types and configurations of LID practices in a mixed land cover watershed using a spatially explicit modeling approach. We simulated 10 scenarios across multiple spatial configurations of LID to evaluate the watershed hydrological responses of three practices—rain gardens (RG), permeable pavements (PP), and riparian buffers (RB)—in a 0.92 km2 watershed with mixed suburban, agricultural, and forest land cover. A spatially-explicit ecohydrological model (VELMA) was used to compare changes in the watershed's water balance before and after LID practice implementation.

Overall, we found that the type and extent of LID practices influence watershed hydrological responses in our study system. Our simulation results indicate that LID practices decreased surface runoff and peak flow, and promoted ET, shallow subsurface runoff and infiltration. However, hydrologic responses and effectiveness varied among LID practices and implementation levels. When LID practices were considered individually, on a LID per unit area basis across all LID implementation levels, RG was more effective in reducing runoff and peak flow, and promoting ET, than PP and RB. However, our results indicated that the 100% implementation of RG was more effective at reducing peak flows during small storms than larger ones, suggesting that LID storage capacities are reduced due to soil saturation during and following large events. Further, both RG and PP increased shallow subsurface runoff and infiltration to almost the same extent at the watershed outlet and the combined LID scenario resulted in the highest performance by increasing shallow subsurface runoff and infiltration, and evapotranspiration by 21% and 15%, respectively, and reductions in peak flow and surface runoff of 8.5% and 8%, respectively.

We conclude that the spatial configurations and extent of LID practices, as well as the model selection and degree of watershed heterogeneity, might be critical for assessing the hydrological responses of watershed-scale LID implementation and must be considered in future research. Further research is needed to apply different LID configurations within mixed land cover watersheds to better understand LID performance and to evaluate the effect of LID practices on nutrient and sediment loads.

**Author Contributions:** H.E.G., C.T.N. and N.H. designed the project and modeling approach; N.H. performed, and R.B.M., B.L.B., A.F.B., K.S.D., J.J.H. and P.P.P. assisted with, the model set up, simulation, and calibration; N.H., H.E.G. and B.P.B. analyzed the results and developed the manuscript's structure; all authors contributed to writing the paper.

**Acknowledgments:** The authors would like to thank Amr Safwat with the CBI Federal Services and Amy Prues with Pegasus, Inc., in Cincinnati, Ohio, and Joong Lee in the Center for Urban Green Infrastructure Engineering, Milford, Ohio, for providing base maps, LID practice configurations in ArcGIS, and impervious area maps.

**Conflicts of Interest:** The authors declare no conflict of interest. Findings and conclusions in this article are those of the authors and do not necessarily reflect or represent the views or policies of the US Environmental Protection Agency. Any mention of trade names, products, or services does not imply an endorsement by the US Government or the US Environmental Protection Agency. The EPA does not endorse any commercial products, services, or enterprises.

#### **References**


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

## *Article* **Modeling Landscape Change Effects on Stream Temperature Using the Soil and Water Assessment Tool**

#### **Mamoon Mustafa 1,\*, Brad Barnhart 2, Meghna Babbar-Sebens <sup>1</sup> and Darren Ficklin <sup>3</sup>**


Received: 13 July 2018; Accepted: 15 August 2018; Published: 27 August 2018

**Abstract:** Stream temperature is one of the most important factors for regulating fish behavior and habitat. Therefore, models that seek to characterize stream temperatures, and predict their changes due to landscape and climatic changes, are extremely important. In this study, we extend a mechanistic stream temperature model within the Soil and Water Assessment Tool (SWAT) by explicitly incorporating radiative flux components to more realistically account for radiative heat exchange. The extended stream temperature model is particularly useful for simulating the impacts of landscape and land use change on stream temperatures using SWAT. The extended model is tested for the Marys River, a western tributary of the Willamette River in Oregon. The results are compared with observed stream temperatures, as well as previous model estimates (without radiative components), for different spatial locations within the Marys River watershed. The results show that the radiative stream temperature model is able to simulate increased stream temperatures in agricultural sub-basins compared with forested sub-basins, reflecting observed data. However, the effect is overestimated, and more noise is generated in the radiative model due to the inclusion of highly variable radiative forcing components. The model works at a daily time step, and further research should investigate modeling at hourly timesteps to further improve the temporal resolution of the model. In addition, other watersheds should be tested to improve and validate the model in different climates, landscapes, and land use regimes.

**Keywords:** stream temperature; SWAT; Marys River watershed

#### **1. Introduction**

Stream temperature is an important water quality parameter that affects physical and chemical processes in streams [1]. Higher stream temperatures in river systems represent a growing concern worldwide and can affect the habitat and life spans of fish [2,3]. According to Eaton and Scheller [4], some fish species will disappear from the water body, if stream temperature transcends an upper limit. Using historical data ranging from 30 to 100 years, Kaushal et al. [5] reported that stream temperatures have been increasing throughout the United States at a rate of 0.009–0.077 ◦C/year, with a significant increase in the western United States. In particular, stream temperatures in the Pacific Northwest have reached historical records—at times, they have exceeded the lethal limit of 21.1 ◦C for some aquatic species such as salmon. For example, in the summer of 2015, the river temperature in the Columbia River in the State of Washington reached the level of 24.5 ◦C and led to the death of 235,000 sockeye salmon out of the total 507,000 that passed through the Bonneville Dam [6].

Similarly, Marys River (Hydrologic Unit Code 17090003) is a tributary of the Willamette River in Northwest Oregon and has experienced increasing temperatures over the last few years. The Oregon Department of Fish and Wildlife (ODFW) conducted a stream survey on flow conditions between the years of 1991 and 1993 and observed that the maximum stream temperature reached the maximum limit of 17.8 ◦C (ODFW). However, recently, increasing levels of human activities have resulted in even higher water temperatures. For example, high water temperatures between 21.1 ◦C and 26.7 ◦C were observed from June to August for some tributaries of Marys River (Marys River Watershed Council), which has resulted in several United States Environmental Protection Agency (USEPA) 303(d) listings for temperature exceedances. These trends may be related to changing climatic drivers as well as land use practices (e.g., harvest of timber, increasing barren lands and clear-cutting areas throughout the watershed) and landscape changes (e.g., urban and agricultural development).

Amidst large-scale landscape and land use changes, preservation of riparian buffers can increase stream shading, thereby helping regulate water temperature along stream reaches [7]. Stream shading intercepts and absorbs a large portion of solar radiation before it reaches the water surface, resulting in less thermal energy that reaches and is stored in streams, which indirectly helps to cool stream temperatures. Brown et al. [8] conducted a study in the Alsea watershed in Oregon along the coast range to study the impact of shading on stream temperatures before and after clear cuts in the watershed. They found that clear-cutting resulted in stream temperature increases of 7.8 ◦C one year after the cuts. Bond et al. [9] investigated the impact of riparian reforestation on summer stream temperatures in the Salmon River in northern California, and they found that partial reforestation lowered stream temperatures by 0.11–0.12 ◦C/km and by 0.26–0.27 ◦C/km for full reforestation.

Since increasing water temperatures have remained a major concern in many watersheds, many models have been proposed to simulate stream temperatures at time scales varying from minutes to months (see Ficklin et al. [10] for a brief review). For this paper, we specifically focus on a semi-distributed mechanistic watershed model called the Soil and Water Assessment Tool (SWAT) [11], which has been extensively used to evaluate the effects of landscape and land use changes on different hydrologic components. Ficklin et al. [10] improved the original stream temperature model within SWAT, which was a linear regression model by Stefan and Preud'homme [12] that correlates thermal energy exchange of air temperature to water temperature. Ficklin et al. [10] developed a daily-scale model for stream temperature prediction by integrating multiple climate and hydrological components, including snowmelt, surface runoff, lateral flow, groundwater flow, and finally air temperature. Several studies have found the Ficklin et al. model [10] produces more realistic simulation results compared to the linear regression proposed by Stefan and Preud'homme [12] (Barnhart et al., 2014; Ficklin et al., 2012; Ficklin et al., 2014 [10,13,14]). However, the model developed by Ficklin et al. [10] does not explicitly account for the different types of radiation that affect thermal energy of water systems, and only recent work has attempted to improve the model by incorporating select radiative components [15]. In general, thermal energy added and removed from any water system consists of incoming radiation that adds thermal energy to the water and results in increasing stream temperatures. This incoming radiation mainly consists of solar radiation coming from the sun, atmospheric longwave radiation, landcover longwave radiation, convection, and evaporation. In contrast, backscattered radiation removes thermal energy and helps to cool temperatures. This radiation consists of emitted longwave radiation from the water surface as well as convection and evaporation. Convection and evaporation radiative components can either add or remove energy from any water body depending on stream temperatures and climate conditions, specifically air temperature, humidity, and wind speed.

This paper highlights model development to explicitly incorporate multiple thermal radiation components into the Ficklin et al. [10] stream temperature model within SWAT. These thermal radiation equations are used within the widely used HEATSOURCE model [16], but until now, these equations have yet to be incorporated into SWAT. HEATSOURCE differs from SWAT because it is a reach-based stream temperature model, whereas SWAT is a watershed model. This means that SWAT simulates upland processes in addition to in-stream processes. HEATSOURCE can model stream temperatures at hourly timesteps and requires site-specific data (e.g., shading, canopy structure, stream morphology) that is oftentimes not available over the entire spatial extent of watersheds. Conversely, SWAT simulates hydrologic components and stream temperatures throughout a watershed using a daily time step and utilizes generally obtainable input data, such as spatially distributed precipitation, temperature, elevation, land use, and soil type. Users may prefer to use SWAT instead of HEATSOURCE when site-specific data is unavailable or when the study goal is to determine the effect of alternative land management scenarios on stream temperatures throughout large, heterogeneous watersheds.

The paper is organized as follows. First, the study area and the SWAT model setup are described. Then, three different SWAT stream temperature models are analyzed, including Stefan and Preud'homm [12] air temperature regression, a mechanistic model by Ficklin et al. [10], and an extension of the Ficklin et al. [10] model in which we specifically incorporate radiative components. We calibrate SWAT for hydrologic discharge in the Marys River watershed, and we compare the three stream temperature models to examine their relative performance for multiple sub-basins with different land use/land cover. We demonstrate the utility of our results by comparing simulations for sub-basins within primarily forested and agricultural landscapes.

#### **2. Methodology**

#### *2.1. Study Area*

The Marys River watershed, shown in Figure 1, is located in the Pacific Northwest of the United States (Hydrologic Unit Code (HUC) 17090003) and is part of the Willamette River basin (HUC 170900) in Oregon. It is one of five major river systems located on the western side of the Willamette River. The area of the watershed is of 782 square kilometers. The highest point of the watershed is at Marys Peak at an elevation of 1280 meters above sea level, and the lowest point is in Corvallis, Oregon, where Marys River drains into the Willamette River at an elevation of 76 meters above sea level. The climate of the watershed in the winter season is mild and wet, with an average winter temperature of 5 ◦C and rainfall during the winter ranges from 1000 mm downstream of the watershed to more than 2500 mm at the highest elevation upstream of the watershed. In general, the watershed tends to be dry, sunny, and warm throughout the summer (Marys River Watershed Council). It has an average summer temperature of 17.5 ◦C. High rainfall intensity results in high stream discharge during winter and spring and mean annual flows of 12–13 m3/s. However, during the summer, flows are generally very low, and discharge sometimes drops below one cubic meter per second. Base flow is a major contributor to the flow of the river, where 61–70% of the total stream flow comes from groundwater contributions [17].

The watershed is divided into three different land use categories: Forest, agricultural, and urban. Most of the watershed (65% of the total area) consists of forest, which is largely located along the western portion of the watershed. In these mixtures of deciduous and evergreen forests, small streams flow over beds of gravel and cobbles with high velocities due to steep slope gradients. Flow leaving the forested region then enters agricultural land in the Willamette Valley that consists mainly of cultivated crops, hay, pasture, wheat, and grass seed production. The streams in this region flow on sand and silt with mild slope gradients, resulting in decreased flow velocities. Furthermore, urban areas are situated further downstream within the Willamette Valley (e.g., the cities of Philomath and Corvallis), and the stream flows over mostly flat to nearly flat gradients. Stream velocities decrease significantly as Marys River enters Philomath, Oregon, and then continues eastwards into Corvallis, Oregon, until it meets the Willamette River at the lowest point located in the watershed.

**Figure 1.** (**a**) Overview of Marys river watershed, and (**b**) land use in Oregon, USA.

#### *2.2. Soil and Water Assessment Tool (SWAT)*

The Soil and Water Assessment Tool (SWAT) is a semi-distributed watershed model that is designed to predict the impact of management on water, sediment, and agricultural chemical yields in gauged and ungauged watersheds [11]. In this study, SWAT was used to simulate the hydrologic and stream temperature dynamics within the Marys River watershed and to evaluate a model extension to the stream temperature model developed by Ficklin et al. [10]. The Marys River watershed was divided into smaller sub-watersheds using pre-defined drainage boundaries and a 10-meter digital elevation model, and then the sub-watersheds in the SWAT model were further divided into smaller units called hydrologic response units (HRUs) using ArcSWAT, a toolbox within ArcGIS for SWAT, with a HRU percentage threshold of 5%. Each HRU is a unique combination of land-use, soil type, and topographic slope and represents the basic unit for conducting mass balances and hydrologic flow in SWAT. The area of the 46 predefined sub-basins varies from 35 square kilometers for the largest sub-basin to 3.0 square kilometers for the smallest. The average sub-basin area is 17 square kilometers. The watershed slope was divided into two categories: (1) A steep gradient area located mostly within the forested regions in the western portion of the watershed, and (2) the nearly flat region located within the Willamette valley, east of the watershed where the cities Philomath and Corvallis are located.

SWAT's input data types include spatial GIS input files such as a Digital Elevation Model (DEM), a land use land cover layer, and a soil layer [18]. Input data needed to delineate the watershed including the DEM, sub-basins, and stream layers in addition to necessary land use and soil SSURGO (Soil Survey Geographic Database) layers to build the HRUs were acquired from United States Department of Agriculture [19]. Three weather stations were used as climate forcings: Corvallis Water Bureau (CWB) COOP ID of (351877), Hyslop weather station, which is also known as Oregon State University weather station (OSU) COOP ID of (351862), and finally Corvallis municipal airport (KCVO) weather station. Weather data of the Corvallis Water Bureau (CWB) and Hyslop weather station were obtained from National Oceanic and Atmospheric Administration (NOAA) for 2005 to 2014. Weather data for Corvallis Water Bureau station included only precipitation and minimum and maximum air

temperature. The Hyslop weather station and Corvallis municipal airport stations included data for precipitation, minimum and maximum temperature, wind speed, and humidity for the period of 2005 to 2014. SWAT was used to simulate flow and stream temperature throughout the Marys River watershed for the period 2005–2014. This includes the period 2010–2014 when observations for stream temperature were available.

#### *2.3. Stream Temperature Models*

#### 2.3.1. Model 1: Linear Regression

The default SWAT stream temperature model uses a linear relation between air temperature and stream temperature developed by Stefan and Preud'homme [12] to calculate stream temperature in the Mississippi River basin, as shown in Equation (1):

$$\mathbf{T\_{water}} = \begin{array}{c} \mathbf{5.0} \ + \ \mathbf{0.75} \times \mathbf{T\_{air}} \end{array} \tag{1}$$

Twater is the average daily water temperature (◦C), and Twater is the average daily water temperature (◦C). Stream temperatures predicted from the above equation will always be higher than air temperature, which is generally a fair assumption for small streams with shallow water depths where stream temperature is primarily controlled by air temperature. However, this may not be necessarily true for streams influenced by snowmelt, surface runoff, and groundwater contributions [10].

#### 2.3.2. Model 2: A Mechanistic Approach Involving Air Temperature and Hydrological Flows

Ficklin et al. [10] developed a mechanistic stream temperature model within SWAT by combining air temperature (heat exchange) and hydrological inputs (flow mixing) including different hydrological parameters, surface runoff, lateral flow, snowmelt, and groundwater contributions. The Ficklin et al. [10] stream temperature model discretizes stream temperature determination into three components: (1) Within the sub-basin, (2) contribution of upstream sub-basins to the targeted sub-basin, and (3) finally heat exchange between air temperature and the stream.

The first part of the stream temperature calculation within the sub-basin (Equation (2)) calculates the local temperature based on a mixing of surface runoff, lateral flow, groundwater, and snowmelt temperatures within the sub-basin flowing to the main stream:

$$T\_{\rm w,local} = \frac{\alpha (0.1 \, \text{Sub}\_{\rm snow}) + \beta \left(T\_{\rm gw} \, \text{Sub}\_{\rm gw}\right) + \lambda \left(T\_{\rm air,lag} \, \text{Sub}\_{\rm surq} + \text{Sub}\_{\rm latq}\right)}{\text{Sub}\_{\rm wyld}}\tag{2}$$

Tair,lag is average daily air temperature with a lag (◦C), and α, β, and λ are calibration coefficients that relate the relative contribution of the hydrologic components to local water temperature (dimensionless). Subsnow is the snowmelt contribution in sub-basin (m3/day), Subgw is the groundwater contribution in sub-basin (m3/day), Subsurq is the surface runoff in the sub-basin (m3/day), Sublatq is the lateral soil flow in sub-basin (m3/day), and Subwyld is the water yield in the sub-basin combining all of the above hydrological inputs (m3/day).

The second part of the Ficklin et al. [10] calculates the effect of upstream sub-basin flow on stream temperature, as shown in Equation (3):

$$T\_{\text{waterinitial}} = \frac{(T\_{\text{w,upstream}})(\text{Q}\_{\text{outlet}} - \text{Sub}\_{\text{wydd}}) + (T\_{\text{w,local}} \times \text{Sub}\_{\text{wydd}})}{\text{Q}\_{\text{outlet}}} \tag{3}$$

Twaterinitial is the stream temperature adding the effects of flow within the sub-basin, Tw,local was calculated previously, Tw,upstream is the water temperature of streams entering the sub-basin (◦C), and Qoutlet is the stream flow discharge at the outlet of sub-basin (m3/day).

*Water* **2018**, *10*, 1143

The final step is to calculate the stream temperature by including the effect of air temperature:

$$\begin{aligned} \text{T}\_{\text{water}} &= \text{T}\_{\text{water} \text{inital}} + \text{K}(\text{T}\_{\text{air}} - \text{T}\_{\text{water} \text{inital}})(\text{TT}) \text{ if } \text{T}\_{\text{air}} > 0\\ \text{T}\_{\text{water}} &= \text{T}\_{\text{water} \text{inital}} + \text{K}((\text{T}\_{\text{air}} + \varepsilon) - \text{T}\_{\text{water} \text{inital}})(\text{TT}) \text{ if } \text{T}\_{\text{air}} < 0 \end{aligned} \tag{4}$$

Twater is the final stream temperature of water (◦C) for a given sub-basin, Tair is the average daily air temperature (◦C), K is the bulk coefficient of heat transfer (1/h), TT is travel time of water through the sub-basin (hour), and finally ε is air temperature addition coefficient (for when air temperature drops below zero).

This mechanistic stream temperature model requires calibration coefficients α, β, γ, k, lag as well as annual groundwater temperatures as inputs. All of the other inputs needed to run the model are provided by SWAT. Groundwater temperature can be estimated from weather data provided as the annual average air temperature, and it is often taken 1–2 ◦C higher than the average annual air temperature [20].

2.3.3. Model 3: A Mechanistic Approach Involving Air Temperature, Hydrological Flows, and Radiative Components

Changes in stream temperature are affected by heat and mass transfers [16] that are dependent on channel morphology, hydrology, and stream vegetation, which provides shading near streams. Vegetation especially helps in cooling temperatures by intercepting and absorbing incoming solar radiation. The mechanistic model introduced in the last section accounts for flow transfer and mixing of various hydrologic components, including sub-basin surface runoff, lateral flow, snowmelt, and ground water, which oftentimes help to reduce temperatures, depending on the season. However, the model utilizes a bulk heat coefficient to account for radiative heat exchange between the air-water interface and does not account for land cover or vegetation near streams. The dependency of the model on the air-water correlation of heat exchange can thereby lead to over-prediction of stream temperatures. Water temperature change related to heat transfer is a function of several sources of radiative heat exchange, as shown in Equation (5):

$$\Phi\_{\text{total}} = \Phi\_{\text{SR}} + \Phi\_{\text{longware-atmosphure}} + \Phi\_{\text{longware-atmosov}} + \Phi\_{\text{convection}} + \Phi\_{\text{conduction}} + \Phi\_{\text{evapation}} \tag{5}$$

where, Φtotal is the net radiation exchange and is equal to the direct and diffuse solar radiation ΦSR as well as the longwave-atmosphere, longwave-landcover, convection, and evaporation components.

Direct and diffusive solar radiation represent the largest sources of incoming thermal energy into streams. Longwave radiation from different sources also plays a role in increasing and decreasing temperatures: Atmospheric and land-cover longwave radiation add energy to water volumes and increase water temperature, while longwave radiation from the water surface emits radiation from the water surface to the atmosphere to cool streams. Energy loss due to evaporation is considered a larger contributor of decreasing stream temperatures when the required energy is met to change the water phase from liquid to gas. Overall, convection or the air-water interface is considered a very small portion of the total energy budget. Groundwater flux also helps to decrease stream temperatures when added as a thermal input.

These radiative components have been included in the HEATSOURCE model but have not been incorporated explicitly within spatially distributed watershed models such as SWAT. Each of the components are defined as follows.

$$
\Phi\_{\rm SR} = H\_{\rm day} - 0.5 \times H\_{\rm day} (1 - e^{-\mathbf{k} \times \rm LAI}) \tag{6}
$$

ΦSR is the amount of solar radiation reaching the water surface, Hday is the incident total solar radiation per day (MJ/m2·day), k is the light extinction coefficient, and LAI is the leaf area index. Atmospheric solar radiation is calculated as follows:

$$
\Phi\_{\text{longware-atmospheric}} = 0.96 \times \varepsilon\_{\text{atm}} \times \sigma \times (\text{T}\_{\text{air}} + 273.15)^4 \tag{7}
$$

Φlongwave-atmosphere is the longwave radiation emitted from the atmosphere, εatm is the emissivity of the atmosphere (unitless), <sup>σ</sup> is the Stefan-Boltzmann constant (MJ·m−2·day−1·K−4), and Tair is the average air temperature per day (◦C). The atmospheric solar radiation depends on the emissivity of the atmosphere (εatm), in which the percent of cloudiness and type of landuse can increase or decrease atmospheric emissivity. The emissivity of the atmosphere is calculated as εatm = 0.767 × (ea) 1 7 , where ea is the vapor pressure of air (mbar) (i.e., H × es), es is the saturation vapor pressure (mbar) (i.e., 6.1275 × e 17.27 <sup>×</sup> Tair 237.3+Tair ), and H is the relative humidity (unitless).

The longwave radiation emitted from landcover Φlongwave-landcover is dependent on the view to the sky θVTS (unitless) and can be calculated as follows:

$$
\Phi\_{\text{longware-landover}} = 0.96 \times \left(1 - \theta\_{\text{VIS}}\right) \times 0.96 \times \sigma \times \left(\text{T}\_{\text{air}} + 273.15\right)^4 \tag{8}
$$

The last component of longwave solar radiation is the water surface longwave:

$$
\Phi\_{\text{longway\textquotedblleft surface\textquotedblright surface}} = \varepsilon\_{\text{w}} \times \sigma \times (\text{T}\_{\text{s}} + 273.15)^4 \tag{9}
$$

Φlongwave-water surface surface is the longwave radiation emitted from the water surface, ε<sup>w</sup> is the water emissivity taken as 0.97, and Ts is the stream temperature (◦C). The evaporation from the water surface is the most effective component in decreasing the thermal energy stored in water:

$$
\Phi\_{\text{evaporation}} \quad = \mathfrak{p} \times L\_{\mathfrak{e}} \times E \tag{10}
$$

ρ is the density of water (kg/m3), Le is the latent heat of vaporization (MJ/kg), which is calculated as Le <sup>=</sup> 2.501 <sup>−</sup> 2.361 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>×</sup> Tair. <sup>E</sup> is the evaporation rate of the water surface (m/day) and is calculated using a mass transfer method: f(w) × (esw − eaw), where f(w) is the wind function <sup>a</sup> + <sup>b</sup> × <sup>w</sup> that depends on coefficients a and b (mbar<sup>−</sup>1) and the wind speed w (m/s) measured 2 m above the water surface [16]. Finally, esw is the saturation vapor pressure of water (mbar), and eaw is the vapor pressure of water (mbar).

The convection radiation component is calculated using the previously calculated evaporative flux Φevaporation and Bowen's ratio BR:

$$
\Phi\_{\text{convection}} = \text{BR} \times \Phi\_{\text{evaporation}} \tag{11}
$$

Here, BR is unitless (i.e., 0.00061 <sup>×</sup> PA <sup>×</sup> Twater <sup>−</sup> Tair (esw <sup>−</sup> eaw) ), and PA is the adiabatic air pressure (mbar) [i.e., 1013 − 0.1055 × z]. z is the measurement height in meters [i.e., >zd + zo]; zd is the zero-plane displacement (m) [i.e., 0.7 × HLc], zo is the roughness height = 0.1 × HLc, and HLc is the height of emergent vegetation (m).

The change of stream temperature due to thermal energy flux is calculated as follows:

$$\mathbf{T\_{W-TD}} = \frac{\Phi\_{\text{total}}}{\rho\_{\text{W}} \times \mathbf{C\_{W}} \times \mathbf{d\_{W}}} \tag{12}$$

Here, Tw−TD( ◦<sup>C</sup> day ) is the temperature change generated from the thermal components, Φtotal is the net driver (MJ/m2·day), Cw (MJ/kg·C) is the specific heat capacity of water, and dw (m) is the depth of water in the channel, which is estimated by the SWAT model.

The final stream temperature is calculated by replacing the second term in Equation (4) by the new generated stream temperature, as follows:

$$\mathbf{T\_{water}} = \mathbf{T\_{water}} \mathbf{i} \mathbf{i} + \mathbf{T\_w} - \mathbf{T} \mathbf{D} \tag{13}$$

Note that, as mentioned above, the majority of the components are calculated within the SWAT model automatically.

Overall, this model is useful because it utilizes all of the distributed information (mechanistically simulated via SWAT) regarding stream temperature and its relation to hydrologic components within a networked watershed, yet it also explicitly incorporates radiative energy exchange at the surface of the stream.

#### *2.4. Model Calibration/Validation Methodology*

SWAT was used to simulate daily hydrologic discharge at each of the sub-basins within the Marys River watershed from 2010–2014 in order to match observed discharge and stream temperature data. The SWAT model was manually calibrated for stream flow between 2010–2014 using the United States Geological Survey (USGS) (14171000) Philomath flow gauge, which is located 6.7 km southwest of where the Marys River meets the Willamette River, covering a 394 km2.

The model was only calibrated without validation due to the availability of observations for flow since the available observations only included a period of less than 10 years. Based on the data availability, the model was calibrated for the period of January 2010 to December 2014.

The Nash Sutcliffe efficiency (NSE; Nash and Sutcliffe (1970)) criterion Equation (14) and Pearson's product moment correlation coefficient (1999; Equation (15)) were used to evaluate hydrologic model efficiency:

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} \left(\text{O}\_{\text{i}} - \text{S}\_{\text{i}}\right)^{2}}{\sum\_{i=1}^{n} \left(\text{O}\_{\text{i}} - \text{O}\_{\text{avg}}\right)^{2}} \tag{14}$$

$$\mathbf{R}^2 = \left(\frac{\sum\_{i=1}^n \left(\mathbf{S\_i} - \mathbf{S\_{avg}}\right) \left(\mathbf{O\_i} - \mathbf{O\_{avg}}\right)}{\left[\sum\_{i=1}^n \left(\mathbf{S\_i} - \mathbf{S\_{avg}}\right)^2\right]^{0.5} \left[\sum\_{i=1}^n \left(\mathbf{O\_i} - \mathbf{O\_{avg}}\right)^2\right]}\right)^2\tag{15}$$

O is the observed value, S is the model prediction, Oavg is the overall observed mean, and Savg is the overall simulated mean. The NSE values range from −∞ to one; a NSE value of less than 0.5 designates an 'unsatisfactory' model, while a NSE value above 0.75 is considered a 'very good' model [21]. R<sup>2</sup> values range from zero to one, with zero indicating a nonlinear relationship between the observed and predicted value and one indicating a perfect fit and a linear relationship between the observed and simulated variables.

The stream temperature models were manually calibrated using root mean square error (RMSE) values as well as percent bias (PBIAS). According to Chai et al. [22], RMSE is widely used as a statistical metric tool to assess performance of models. RMSE can be calculated as follows:

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} \left(\text{O}\_{i} - \text{S}\_{i}\right)^{2}}{n}} \tag{16}$$

Where O is observed value, S is model predicted value, and n is the total number of the points. A RMSE value of zero indicates a perfect fit. According to Singh et al. [23], values of RMSE less than half of the standard deviation of the measured data can be taken as acceptable for the model evaluation.

In addition to RMSE, PBIAS was also used to evaluate the models. PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts [24]. Zero is the optimal value of PBIAS, and a low absolute value implies an accurate model. According to Gupta et al. [24], positive PBIAS values indicate a model underestimation bias, and negative PBIAS values indicate a model overestimation bias. PBIAS can be calculated as follows:

$$\text{PBIAS} = \frac{\sum\_{i=1}^{n} \left(\text{Oi} - \text{S}\_{\text{i}}\right) \times 100}{\sum\_{i=1}^{n} \left(\text{Oi}\right)} \tag{17}$$

A set of seven calibration parameters were selected and manually modified to calibrate SWAT for hydrology (Table 1), and five parameters were manually modified to calibrate SWAT for stream temperature (Table 2). Observed daily stream temperature data corresponding to SWAT's sub-basins 8, 15, and 17 were available from 2010 through 2014, while observations for sub-basin 36 extended from 2011 through 2014.


**Table 1.** Streamflow calibration parameters.

\* Values are percentages of the original values.


To compare the differences in model performance, kernel density estimates were calculated using R software. This nonparametric technique is similar to using histograms to highlight the differences between model simulations and observed data for the three tested models.

#### **3. Results and Discussion**

#### *3.1. Hydrology Calibration in SWAT*

As mentioned previously, SWAT was used to simulate daily hydrologic discharge at each of the sub-basins within the Marys River watershed for 2005–2014, which included 2010–2014, when observations for stream temperature were also available. The NSE and R<sup>2</sup> values of the default model's simulated flow (no calibration) were −0.37 and 0.50, respectively, indicating an unsatisfactory model. From Figure 2, it is clear that the uncalibrated model overpredicts the peaks during storm events. Also, the model's responses to each rain event are very rapid, and the water loss rates are excessive, which results in zero flow for late summer periods.

**Figure 2.** (**a**) Uncalibrated daily SWAT discharge simulations compared with observations (2010–2014) from the USGS (14171000) Philomath flow gage. (**b**) Uses a log scale to emphasize low-flow conditions.

Manual calibration resulted in model improvement of NSE values from −0.37 to 0.72. This designates a 'good' model according to Moriasi et al. [21], since it is >0.65. The R2 value of the model increased to 0.80. Figure 3 shows that the calibrated model matches both the base flow and peaks well, whereas all of the peaks were mainly over predicted by the default model. The default model could not capture the very low stream flows for some summer days and resulted in zero flow (Figure 2), but the calibrated model (Figure 3) fixed the low-flow problem and improved the results. Figure 2 shows that the maximum simulated peak for the uncalibrated model was around 500 m3/s, whereas the calibrated model (Figure 3) reduced this value to match the peak observed around 300 m3/s. The rapid response of the main channel to any storm event led to the over-prediction of peaks even for small rain storms in the uncalibrated model. Also, the lateral flow travel time parameter helped to slow the response to storm events and smoothed the hydrograph. The other major problem was associated with the excessive loss of water in a short period of time after a rapid response to any storm; water in the uncalibrated model was lost instantaneously and therefore resulted in zero flow for late summer days. The effective hydraulic conductivity in the main channel (CH\_K2) was used to prevent the excessive loss from the main channel and managed to eliminate the zero flow days. Also, this parameter helped to smooth and eliminate the transient fluctuations in the hydrograph. The SCS curve number for moisture conditions II (CN2) was used to reduce or increase the simulated peaks to match the observed hydrograph. The canopy interception parameter for specified land cover types (CANMX) as well as the minimum and maximum snowmelt factors (SMFMN, SMFMX) all helped to increase the evaporation rate. The high surface runoff surge was decreased using the soil evaporation compensation factor

(ESCO), the SCS curve number for moisture conditions II (CN2), and the threshold depth of water in the shallow aquifer (GWQMN).

**Figure 3.** (**a**) Manually calibrated daily Soil and Watershed Assessment Tool (SWAT) discharge simulations compared with observations (2010–2014) from the United State Geological Survey (USGS) (14171000) Philomath flow gage. (**b**) Uses a log scale to emphasize low-flow conditions.

#### *3.2. Stream Temperature Calibration in SWAT*

After a satisfactory hydrologic calibration was performed, manual calibration was performed for two of the three stream temperature models that will be tested in the Marys River using SWAT. The first model is the linear regression model from Stefan and Preud'homme [12] and was not calibrated. The second model is the Ficklin et al. [10] model, and the third model is our extension to the Ficklin et al. [10] model. As shown in Table 2, we manually calibrated five parameters to best match the second and third models to the stream temperature observations between 2010 and 2014. Note that only summer stream temperature observations were available and that only a single set of calibration parameters were used for the entire watershed.

#### *3.3. Stream Temperature Model Comparison*

After satisfactory hydrologic and stream temperature calibrations were performed, three models were tested to simulate stream temperature in the Marys River using SWAT. The first was the linear regression model from Stefan and Preud'homme [12], the second the Ficklin et al. [10] model, and the third is our extension to the Ficklin et al. [10] model that incorporates HEATSOURCE radiative forcing components to the Ficklin et al. [10] model. For the remainder of this paper, these models will be referred to as Model 1, Model 2, and Model 3, respectively.

Table 3 shows RMSE and PBIAS results for these three models using daily data. Kernel density estimates for the differences between simulated and observed stream temperatures for the four sub-basins are plotted in Figure 4.


**Table 3.** Comparison of root mean square error (RMSE) and percent bias (PBIAS) values for the three tested stream temperature models (daily).

**Figure 4.** Kernel density estimates for the differences between observed and simulated stream temperatures as calculated using the three models for four sub-basins (**a**–**d**). The vertical line at zero indicates simulations that exactly match the observed data.

In general, Model 1 performed the worst among all three models (see Table 3 and Figure 4), and Model 2 outperformed Models 1 and 3 for all sub-basins. Model 1 consistently overestimated stream temperatures for all of the sub-basins. This is apparent in Figure 4, where the distributions of the differences between the simulated and observed temperatures are shifted from zero for Model 1. The coefficients of Model 1 guarantee that the simulated stream temperature will be above average daily air temperature values when the average air temperature is less than 20 ◦C. Therefore, inclusion of cold water from groundwater or upstream sources is not captured in this model, thus resulting in overestimations. Model 2, which is the calibrated Ficklin et al. [10] model, showed improvements compared to the linear model (Model 1) in both Table 3 and Figure 4, which agrees with previous studies (Barnhart et al., 2014; Ficklin et al., 2012; Ficklin et al., 2014) [10,13,14]. This is presumably because Model 2 incorporates hydrologic components, including groundwater upstream temperatures, in addition to an air-heat exchange transfer coefficient. Model 3, which replaced the simple air-heat exchange transfer coefficient from Model 2 with explicitly calculated radiative components, shows similar distributions with Model 2 (Figure 4), yet the performance values of Model 2 are

generally better. This is likely due to the high variability of the radiative components included in Model 3, which will be discussed further in the next section.

#### *3.4. Land Cover Effects on Stream Temperature*

We now compare Models 2 and 3 to demonstrate that the incorporation of radiative components is able to simulate the influence of land use and land cover on stream temperatures (e.g., forested vs. agricultural regions). Stream temperature simulations for two sub-basins—a forested area (sub-basin 8) and an area dominated by agriculture with low vegetation (sub-basin 36)—are shown in Figure 5.

**Figure 5.** Stream temperatures simulated using the Ficklin et al. [10] model (Model 2; panel (**a**)) and the extended version that includes radiative components (Model 3; panel (**b**)). Model 2 does not capture changes in stream temperature associated with land cover type, while Model 3 simulates increased temperatures for agricultural regions.

Model 2 simulates nearly identical stream temperatures for both forested and agricultural sub-basins. Conversely, by including the various radiative components into the model, Model 3 simulates consistently increased stream temperatures associated with the agricultural sub-basin. This is especially apparent during the early summer months. To examine differences in the models further, Figure 6 compares the net radiative components of Model 3 Equation (5) with the K-component second term in Equation (4) of Model 2.

Figure 6 shows that the net radiation as calculated using Model 3 (orange lines) changes according to the primary land use cover for a given sub-basin in SWAT. For example, agricultural sub-basins have larger incoming (positive) radiation contributions that help to increase stream temperatures, whereas forested areas have small or negative radiative effects, depending on the season, due to increased LAI, reduced solar radiation reaching the water surface, and evaporative fluxes. Model 2 uses a bulk coefficient of heat transfer in the second term of Equation (4). This reflects a convection component of the net radiative balance, but it does account for the other radiative energy terms, including solar radiation, atmospheric longwave, land surface longwave, water surface longwave, and evaporation. Therefore, it is not able to capture cover-related differences in net radiation and therefore changes in stream temperature due to landscape and land use changes.

Figure 6 also shows that the net radiative driver as calculated in Model 3 has much higher variability than the K-component used in Model 2. We found that the high variability (i.e., noise) is mainly due to SWAT's estimation of solar radiation as well as the evaporation calculations, which are not explicitly accounted for in Model 2. This likely led to the reduced model performance exhibited when comparing model simulations to observed data.

**Figure 6.** Comparison of the thermal energy impacts on stream temperatures for Models 2 and 3 in both primarily forested and agricultural sub-basins.

Figure 7 compares these data further by plotting the observed stream temperature as well as simulations using Models 2 and 3 for both the forested (sub-basin 8) and agricultural (sub-basin 36) sub-basins.

**Figure 7.** Comparison of observed and simulated stream temperatures for agricultural and forested sub-basins. The observed data (black) show a consistent increase in stream temperatures for the agricultural sub-basin (sub-basin 36). Model 2 does not capture consistent changes in stream temperature associated with land cover type. Model 3 simulates increased temperatures for the agricultural sub-basin, but it overestimates the effect compared to the observed data.

The mean (standard deviation) differences between the stream temperature simulations for the forested sub-basin (8) and the agricultural sub-basin (36) were −0.68 ◦C (1.30 ◦C) and −4.1 ◦C (2.34 ◦C) for Models 2 and 3, respectively. These simulations can be compared with the difference between the observed data for the forested and agricultural sub-basins (8 and 36): −1.07 ◦C (0.63 ◦C). Note that Model 2 gives estimates that are closer to the observed values than Model 3; this is not surprising since the RMSE values for Model 2 were lower than the values for Model 3 (Table 3). However, Model 2 is symmetric about the 1:1 line shown in Figure 7 and does not capture the positive bias of agricultural stream temperatures that is shown by the observed values as well as in the simulated values of Model 3. Model 3 captures the increases in stream temperature associated with the lack of forest cover in the agricultural sub-basin, yet Model 3 overpredicts this effect, and the variation in the simulated values are much greater than Model 2 or the observed data.

Overall, Model 3 may be useful for simulating the watershed-scale impacts of land use conversion from forest to agriculture on stream temperatures; however, model improvement—potentially through improved calibration—is needed to better match observed data. In addition, while Model 3 is able to simulate impacts of land use on stream temperature, this advantage is produced with the trade-off that radiative components exhibit much higher variability, and this noisy fluctuation (which can be seen directly from Figures 5–7) further decreases RMSE values (Table 3).

#### **4. Conclusions**

This study sought to explicitly incorporate radiative forcing components into an existing mechanistic, semi-distributed stream temperature model Ficklin et al. [10] using SWAT. Ultimately, stream temperature is controlled by different climate components including humidity, wind speed, evaporation, and solar radiation besides air temperature and is also heavily dependent on hydrologic components including surface flow, groundwater, and snow melt processes. Our new model leverages the Ficklin et al. [10] stream temperature model, which accounts for hydrological and climatological components and discretizes the prediction of stream temperature into three parts: (1) Streamflow within a sub-basin, (2) contributing upstream sub-basin hydrologic components, and (3) accounting for the heat exchange between air and water surface, which can be described as a convection term. The extended model replaces the convection K-component term within the Ficklin et al. [10] model with a more comprehensive characterization of radiative energy terms, including solar radiation, atmospheric longwave, land surface longwave, water surface longwave, evaporation, and finally convection drivers. The extended model was used along with the Ficklin et al. [10] model and the linear regression model to simulate stream temperatures within agricultural, forested, and mixed sub-basins within the Marys River watershed. Results showed that all models performed reasonably well, and the Ficklin et al. [10] model outperformed the others. However, the extended model was capable of simulating differences between stream temperatures associated with agricultural and forested watersheds that reflected observed data, although the differences were overestimated. The reduced performance of the extended model that included radiative components might be able to be improved by further calibration; yet, the high variability of the radiative terms is also limiting. For example, the model relies on incoming solar radiation as well as wind velocity, which are difficult to represent over large spatial scales and feature high variability. Alternative formulations of the radiative components should be considered in future work. Bogan et al. [1] suggests using a shading factor instead of a leaf area index, which may lead to more accurate results, assuming shading information is available for the watershed of interest. In addition, alternative estimations for solar radiation, or perhaps any of the radiative energy terms, could improve the model and should be pursued. Overall, incorporating radiative components into the Ficklin et al. [10] stream temperature provides a new mechanism for simulating the effects of alternative land uses on stream temperature within SWAT. This will be especially useful for land managers and decision makers when considering alternative land management scenarios and conservation strategies using SWAT.

**Author Contributions:** Conceptualization, M.B.-S.; Methodology, M.B.-S.; Software, M.M. and B.B.; Validation, M.M. and B.B.; Formal Analysis, M.M.; Writing—Original Draft Preparation, M.M. and B.B.; Writing—Review & Editing, M.B.-S. and D.F.; Visualization, M.M. and B.B.; Supervision, M.B.-S.

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

**Acknowledgments:** The authors thank reviewers for their comments that helped in preparing a better version for publication. Also, we thank the editorial team of Water for proofreading and preparing the paper for publication. The views expressed in this article are those of the authors and do not necessarily represent the views or policies of the U.S. Environmental Protection Agency. Any mention of trade names, products, or services does not imply an endorsement by the U.S. Government or the U.S. Environmental Protection Agency. The EPA does not endorse any commercial products, services, or enterprises.

**Conflicts of Interest:** We declare no conflicts of interest.

#### **References**


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

## *Article* **Improved Soil Temperature Modeling Using Spatially Explicit Solar Energy Drivers**

**Jonathan J. Halama 1,2,\*, Bradley L. Barnhart 1, Robert E. Kennedy 2, Robert B. McKane 1, James J. Graham 3, Paul P. Pettus 1, Allen F. Brookes 1, Kevin S. Djang <sup>4</sup> and Ronald S. Waschmann <sup>1</sup>**


Received: 31 August 2018; Accepted: 4 October 2018; Published: 9 October 2018

**Abstract:** Modeling the spatial and temporal dynamics of soil temperature is deterministically complex due to the wide variability of several influential environmental variables, including soil column composition, soil moisture, air temperature, and solar energy. Landscape incident solar radiation is a significant environmental driver that affects both air temperature and ground-level soil energy loading; therefore, inclusion of solar energy is important for generating accurate representations of soil temperature. We used the U.S. Environmental Protection Agency's Oregon Crest-to-Coast (O'CCMoN) Environmental Monitoring Transect dataset to develop and test the inclusion of ground-level solar energy driver data within an existing soil temperature model currently utilized within an ecohydrology model called Visualizing Ecosystem Land Management Assessments (VELMA). The O'CCMoN site data elucidate how localized ground-level solar energy between open and forested landscapes greatly influence the resulting soil temperature. We demonstrate how the inclusion of local ground-level solar energy significantly improves the ability to deterministically model soil temperature at two depths. These results suggest that landscape and watershed-scale models should incorporate spatially distributed solar energy to improve spatial and temporal simulations of soil temperature.

**Keywords:** soil temperature; solar energy; watershed model; landscape scale; VELMA

#### **1. Introduction**

Soil temperature affects several key ecosystem properties. Through surface runoff and subsurface groundwater transport, soil temperatures can lead to increased stream temperatures, which in turn impact salmonid and other fish habitats [1]. Soil temperatures mediate rates of biogeochemical transformations in soils, strongly influencing local to global-scale patterns in the cycling, retention and loss of carbon and nutrients from ecosystems [2]. Seasonal soil temperature trends can shift photosynthetic recovery timing and therefore impact overall net primary production (NPP) [3]. Such soil temperature effects are subject to modification by physical landscape factors, such as object shading, slope aspect and thermal isolation from a detritus layer or snow pack [4].

Mechanistic watershed models such as Visualizing Ecosystem Land Management Assessments (VELMA) (2.0, U.S. Environmental Protection Agency—Western Ecology Division, Corvallis, OR, USA) [5], Soil and Water Assessment Tool (SWAT) (2009, Texas A&M, College Station, TX, USA) [6], Regional Hydro-Ecologic Simulation System (RHESSys) (2004, University of California, Santa Barbara, CA, USA) [7], and Hydrologic Simulation Program—Fortran (HSPF) (11, U.S. Environmental Protection Agency—National Exposure Research Laboratory, Athens, GA, USA) [8] use a mechanistic (as opposed to statistical) approach to model hydrodynamics throughout a watershed using sub-daily or daily time steps. Models such as these utilize equations to simulate hydrologic dynamics and soil moisture by tracking the rate of water transfer based on soil porosity, soil depth and the available precipitation.

Watershed models simulating soil temperature at multiple depths rely on several observed and simulated environmental variables including air temperature, precipitation, soil moisture, soil depth, and physical soil properties. For national and regional-scale modeling purposes, climate data can often be obtained from various governmental agencies such as the National Oceanic and Atmospheric Administration (NOAA), Natural Resources Conservation Service (NRCS) and others that maintain large-scale networks of climate monitoring stations such SNOTEL and SCAN [9]. At more local scales, climate data collection tends to focus on site-specific requirements, such as municipal airports, Long Term Ecological Research Stations and university research forests [10–13]. This can result in data limitations for spatially explicit models [14]. Several groups process the site data to produce spatial datasets at various spatial scales (e.g., Parameter-elevation Regressions on Independent Slopes Model (PRISM) and Daily Surface Weather Data (Daymet)) [15,16]. Soil column properties can be acquired through field work or obtained by utilizing soil datasets (e.g., State Soil Geographic (STATSGO) and Soil Survey Geographic (SURGO)) [17].

While climate and soil properties influence soil temperatures, solar energy is the most significant environmental variable influencing soil temperature. There are existing spatial models that account for solar energy inputs at a local or stream reach scale: SHADE2 (1.0, University of Georgia, Athens, GA, USA) [18], HeatSource (8.0, Oregon Department of Environmental Quality, Portland, OR, USA) [19], and iLand (1.0, Seidl and Rammer, Vienna, Austria) [20], utilize small-area representations of solar energy. However, these model's spatial heterogeneity of shade is utilized for quantifying shade along stream reaches or within forest plots; not ground-level irradiance within models representing complete watersheds or landscapes.

Previous methods of incorporating solar energy within model representations of complete watersheds employ one or more proxy variables (e.g., canopy coverage or air temperature), or they simply utilize an average daily global irradiance value. Current environmental mechanistic models (e.g., VELMA, SWAT, RHESSys, HSPF) [5–8] use a global solar irradiance subroutine to calculate the total solar energy (W/m2) for the entire watershed or for sub-catchments. Due to their simpler ground-level solar energy representations, although these models capture the seasonal pattern of irradiance, they lack a spatially heterogenous representation of topographic and landscape object shading that affects ground-level solar energy levels. Solar energy estimate methods themselves include uncertainty due to environmental variables like cloud fraction and albedo [21]. Additional variables (i.e., aerosol optical properties, cloud asymmetry, water vapor distribution) may have seasonal and regional influence on the accuracy of solar energy estimates [22].

There remains a gap in soil temperature modeling where current approaches utilize global solar energy models that are not capturing local energy interactions. Finer-scale spatially distributed estimates of ground-level shade or solar energy could provide improved soil temperature estimates. This paper addresses this gap by incorporating local solar energy data within the soil temperature subroutine of an ecohydrological watershed model to determine its effect on simulated soil temperature predictions at multiple depths.

To demonstrate the utility of linking spatially explicit, ground-level solar energy data with an environmental model, we focus on improving VELMA's soil temperature subroutine by incorporating spatially heterogeneous solar energy as an input driver. First, we discuss a commonly used global solar energy model and soil temperature model that are found in many deterministic watershed models, including VELMA. Next, we present VELMA's original soil temperature subroutine that incorporates spatially explicit inputs of soil moisture and air temperature but does not employ any form of solar energy presentation. Following the original model, we then present VELMA's modified soil temperature subroutine that incorporates spatially explicit inputs of ground-level solar energy, along with the inclusion of soil moisture and air temperature. We compare the predictive skill of the original and new forms of VELMA's soil temperature subroutines using multiple United States Environmental Protection Agency (EPA) Oregon Crest-to-Coast Environmental Monitoring Transect (O'CCMoN) sites. This transect consists of several paired sites of forested and open landscapes [23]. We present results that demonstrate the benefit of including spatially explicit representations of solar energy within watershed-scale models that simulate soil temperature.

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

Watershed models typically include solar energy directly or through a proxy variable to facilitate energy requirements needed within subroutine routines. Plant growth models may require a daily input of solar energy reaching the canopy to drive photosynthesis [24]. Stream temperature models predict shifts in water temperature through variables representing landscape shading, water temperatures, and air temperatures, all of which are solar energy proxies. Snowmelt models may need a daily input of solar energy or air temperature to drive snow melt [25]. While all these subroutines rely on solar energy at the earth's surface, mechanistic models generally lack the ability to capture the spatiotemporal dynamics of solar energy reaching the ground post shadowing.

The soil temperature subroutine from VELMA [5] was chosen for testing. VELMA is a spatially distributed watershed model that simulates hydrologic and biogeochemistry processes within a gridded framework under mechanistic cell interactions. Using a gridded framework, VELMA describes each grid cell as having a ground-level surface and four sub-surface voxels representing the landscapes soil strata. Each subsurface voxel is characterized by soil porosity and soil depth. Water transfers at a daily time step through VELMA's voxel framework. Based on water transmission, nutrients and thermal energy migrate through the simulated soil substrate under mechanistic rules.

#### *2.1. Previous Solar and Generalized Soil Modeling Methods*

Watershed models often employ a clear-sky solar energy model when direct solar energy units are required. A common approach, and the method used by VELMA, is to calculate the clear-sky solar energy that reaches the earth's surface as in Equation (1), where *R* is solar irradiance (W/m2), *ecc* is the eccentricity correction factor, *w* is the Earth's constant angular velocity, *T* is the time frequency, *dec* is the solar declination, and *γ* is latitude [26]:

$$R = \left(24/\pi\right) \times 4.921 \times \text{acc} \times \left[wT \times \sin(d\mathbf{c}) \times \sin(\gamma) + \cos(d\mathbf{c}) \times \sin(wT) \times \cos(\gamma)\right] \tag{1}$$

The VELMA model uses Equation (1) to describe the amount of solar energy reaching the troposphere under clear sky conditions. This approach does not account for the topographic or object shading that locally reduces ground-level solar energy.

Soil temperature modeling within many watershed models, whether utilizing a gridded representation of the landscape or aggregating to sub-catchment scales, typically uses some version of the Carslaw and Jaeger equation to quantify seasonal variation in soil temperature [27]:

$$T\_{\rm soll}(z\_{\prime}, d\_{\rm H}) = T\_{AA} + A\_{\rm surf} \times \text{e}^{-z/dd} \times \sin(\omega \times d\_{\rm h} - z/dd) \tag{2}$$

where *Tsoil*(*z*, *dn*) is the soil temperature (◦C) at depth *z* (mm) for day of the year *dn*, *TAA* is the average annual soil temperature, *Asurf* is the amplitude of the surface fluctuations, *dd* is the damping depth (mm), and ω is the angular frequency of the damping oscillations by day (*dn*). At *z* = 0, the soil temperature reduces to the following:

$$T\_{soil}(0, d\_n) = T\_{AA} + A\_{surf} \times \sin(\omega \times d\_n) \tag{3}$$

which is the average soil temperature perturbed by surface temperature fluctuations and reflects seasonal solar patterns. Conversely, at infinite depth, the soil temperature becomes equal to the annual average soil temperature. This formulation provides a relatively simple method for calculating soil temperatures at multiple depths throughout a watershed. However, the model requires specification of soil heat capacity as well as thermal conductivity to correctly specify the amplitude coefficient and the damping depth.

#### *2.2. Soil Temperature Variations Due to Landscape Coverage*

Temperature profiles of soils can dramatically vary between a forested versus open environment, even if the sites are located proximally near one another. Two sites can be exposed to very similar climate conditions, though due to forest canopy shading, the forested site will have reduced air temperature and a reduction in solar energy loading upon the soil surface. Figure 1 shows observed 2005 daily soil temperature differences between open (clear-cut harvested) site data minus forested site data for the O'CCMoN Soapgrass field site in Oregon (see Section 2.5.1 for field site locations).

**Figure 1.** Daily differences in soil temperature (◦C) between open (clear-cut) and forest sites at the O'CCMon Soapgrass site. Temperature differences were calculated as the Open Site temperature minus the forest site temperature at each Julian day during the year 2005. Thus, positive values denote days where the open site soil temperature was warmer than the forest site soil temperature. Open minus forest soil temperature differences were calculated for each of the two soil depths, 15 cm (red line) and 30 cm (blue line) below the soil surface. The data gap between Julian days 67 and 74 was due to sensor errors.

Figure 1 highlights the soil temperature differences between open and forested sites. A positive temperature means the open site temperature was warmer than the forest site; conversely, negative temperature means the open site temperature was colder than the forest site. Two main observations should be made here: (1) layer 1 is always warmer than layer 2, and (2) for both soil layers the open site is always warmer in summer than the forest site, but is comparatively colder in winter and especially so in layer 2. The open site is significantly warmer from Julian day 45 through 310 (14 February through 6 November), with a peak difference of 4.2 ◦C on Julian day 111 (21 April).

Seasonal differences in the warming and cooling of soils in the open and forest sites (Figure 1) certainly reflect changes in air temperature along with some complicating effects associated with inter-site variations in snow pack and associated insulative properties (Figure 2). Other factors undoubtedly also come into play, such as the effects of seasonal changes in soil moisture (dry summers, wet winters) on soil thermal transmissivity.

Nevertheless, observed increases in summer air temperatures for the open site tended to be 2.12 ◦C warmer than the forest site (Figure 2). The lowest thermal difference was −7.65 ◦C while the highest thermal difference was 8.0 ◦C. During the same period, soil layer 1 open site averaged 1.71 ◦C warmer than the forest site with a minimum of 1.0 ◦C and a maximum of 2.6 ◦C. However, soil layer 2 open site averaged −0.16 ◦C cooler than the forest site with a minimum of −1.0 ◦C and maximum of 1.3 ◦C. That is, while air temperature differences are driven by differences in solar radiation, in the open site there is clearly an additional direct effect of solar radiation on heat transfer to the ground surface and consequent warming of the soil column. This observation underscores the importance of quantifying the direct effect of solar radiation on soil temperature, in combination with effects of soil moisture and other factors mentioned above.

**Figure 2.** Daily air temperature (◦C) and snow depth (cm) for both the forest and open sites at the Soapgrass station. A positive air temperature difference means the open site was warmer than the forest site. A positive depth difference means the open site had more snow than the forest site.

#### *2.3. Original VELMA Air Soil Temperature (AST) Subroutine*

For calculating spatially distributed soil temperature, VELMA accounts for soil moisture damping and the oscillatory effects of solar energy through a modified version of the Carslaw and Jaeger equation. This approach accounts for the seasonal solar energy variability through a time phase lag modification of observed air temperature combined with a temperature modification based on a soil depth attenuation. VELMA's subroutine, along with the equation previously presented by Carslaw and Jaeger (1959) (Equations (2) and (3)), does not account for spatial heterogeneity of solar energy reaching the ground due to topographic or object shading; instead, VELMA's input variables account for the shift in soil temperature due to daily air temperature, soil moisture and soil depth (Figure 3).

**Figure 3.** Combined AST and AST-Solar subroutines' schematic depicts the input data driving both models and displays the crucial data input difference between AST and AST-Solar with Temporal Phase Lag influencing only AST and Solar Energy influencing only AST-Solar.

For each layer, the AST subroutine calculates the soil temperature based on the thermal attenuation of daily air temperature. The input variables for the AST and AST-Solar models are mostly the same but might be utilized differently within each subroutine's equation setup, with the exceptions of the AirLAG and ReducerSOLAR (Table 1).


**Table 1.** AST and AST-Solar subroutine input variables.

The degree of attenuation is adjusted daily by the depth and soil moisture of each soil layer. GTEMP is the resulting soil temperature due to: AirAVETEMP being the daily average air temperature (Table 1), AirLAG (Equation (5) from a prior AirAVETEMP (Table 1) based on seasonal oscillation, SoilDAMPING (Equation (7) influencing the soil moisture based on seasonal oscillation, DepthATTENUATION (Equation (9) based on soil depth and soil damping (Equation (7), PhaseLAG (Equation (6) influenced by the AirAVETEMP (Table 1) based on seasonal oscillation driven by LSDEPTH (depth to surface) (Table 1), SoilDAMPING (Equation (7), LTD (Equation (8) being the soil temperature accumulation at depth, and SoilBELOW (Table 1) influencing soil temperature from the lower soil layer:

GTEMP = AirAVETEMP + (AirLAG − AirAVETEMP − SoilDAMPING) × DepthATTENUATION (4)

AirLAG = Past Air Temperature at Julian Day's PhaseLAG (5)

$$\text{Phase}\_{\text{LAG}} = \left(\text{LS}\_{\text{DEPTH}} / \text{Soil}\_{\text{DAMIPNG}}\right) \times \left(\text{365} / 2\pi\right) \tag{6}$$

$$\text{Soil}\_{\text{DAMPING}} = \sqrt{\frac{\text{LTD} \times 365}{\pi}}\tag{7}$$

$$\text{LTD} = \text{LTD}\_{\text{ACUMULATION}} / \text{Soil}\_{\text{BELLOW}} \tag{8}$$

$$\text{Depth}\_{\text{ATTERUATION}} = \mathbf{e}^\* (-\text{LSDEP} \text{TH} / \text{SoilDAMPING}) \tag{9}$$

Each day, the air temperature is given as an input for each cell within VELMA's watershed framework. The SoilDAMPING (Equation (7) and LSDEPTH (Table 1) variables are used to dampen the variations of soil temperature at larger depths. Any soil temperature shifts due to solar energy are incorporated via a proxy of two oscillatory equations driven by past air temperature (Equation (5) and soil moisture damping (Equation (7). VELMA AST subroutine's performances for open and forested landscapes are tested in the model testing section below.

#### *2.4. New VELMA Air Soil Temperature-Solar (AST-Solar) Subroutine*

The previously described soil temperature model does not utilize spatially explicit solar energy data, so we improved the model by adding the capacity to utilize spatially distributed ground-level solar energy to the current VELMA AST subroutine. This new model is called Air Soil Temperature-Solar (AST-Solar). The inclusion of solar energy within VELMA's original AST model was mainly accomplished through the addition of the new parameter ReducerSOLAR (Equation (12). Like a natural system, solar energy, via the variable ReducerSOLAR (Equation (12), only impacts the top soil layer, called NetSoilTEMP1 (Equation (10). NetSoilTEMP1 is defined as the following:

$$\text{NetSoil}\_{\text{TEMP1}} = \text{Air}\_{\text{TEMP}} \times \text{Reducer}\_{\text{SOLAR}} \times \text{Damping}\_{\text{SOL}} \tag{10}$$

where daily average air temperature (AirTEMP; Table 1), DampingSOIL (Equation (11), and ReducerSOLAR (Equation (12) are multiplied together. AST soil moisture damping was included, yet simplified to only the inversion of each layer's fraction of volume to volume (*v*/*v*) soil moisture (LayerSM; Table 1):

$$\text{Damping}\_{\text{SOIL}} = 1 - \text{Layer}\_{\text{SM}} \tag{11}$$

Solar energy was built into the AST-Solar approach by accounting for the proportional relationship between each cell's solar energy to the landscape cell with the maximum solar energy. At each simulation timestep, the spatially distributed solar energy and the watershed's maximum solar energy at any location are both used to calculate the solar energy reduction called ReducerSOLAR (Equation (12). ReducerSOLAR represents the localized reduction in soil temperature due to shadowing in relation to the watershed's maximum ground-level solar energy. The reduction of solar energy is calculated as follows:

$$\text{Reducer}\_{\text{SOLAR}} = 1 - \alpha \times (1 - (\text{Cell}\_{\text{SOLAR}} / \text{Max}\_{\text{SOLAR}})) \tag{12}$$

where CellSOLAR is each cell of interest within the VELMA framework, MaxSOLAR is the landscape's maximum solar energy value amongst all landscape cells per time step, and α is a calibration factor. The calibration factor α is a fraction [0.0–1.0, where 0.0 is no solar energy and 1.0 is no change to the solar energy] that allows control over the influence of ReducerSOLAR. But, to allow a direct and fair comparison of AST to AST-Solar, calibration factor α was not used in these tests.

VELMA utilizes four soil layers, and the thickness of each layer is customizable. VELMA soil temperature is not directly affected by solar energy, but rather through soil depth attenuation. For layers 2, 3 and 4, the soil temperature (NetSoilTEMPX) is calculated using the 2-day running average temperature of the soil layer directly above, plus a reduction by DampingSOIL (Equation (11):

$$\text{NetSoil}\_{\text{TEMPX}} = \text{Soil}\_{\text{AE\\_TEMP}} \times \text{Damping}\_{\text{SOL}} \tag{13}$$

$$\text{Soil}\_{\text{AVE\\_TEMP}} = \text{(SoilTermp}\_{\text{(IDAy)}} + \text{SoilTermp}\_{\text{(Day}-1)})/2 \tag{14}$$

where SoilTemp(JDay) (Table 1) is the current time steps soil temperature, and SoilTemp(JDay−1) (Table 1) is the prior time steps soil temperature. The soil moisture is applied as the DampingSOIL coefficient (Equation (11).

#### *2.5. Subroutine Testing*

We utilized data from EPA's O'CCMoN sites to test any change in accuracy and seasonal performance between the AST versus AST-Solar soil temperature subroutines. The O'CCMoN transect dataset provided observed driver data of air temperature, photosynthetic active radiation (PAR) as micromoles/meter2/second (μmoles/m2/s), and soil moisture as volume to volume at two soil layer depths [23]. EPA's observed O'CCMoN data helped to compare the AST versus AST-Solar models. Each O'CCMoN site also provided observed soil temperature at two depths. The soil temperature data were used to generate goodness-of-fit metrics against the simulated model results.

#### 2.5.1. O'CCMoN Testing Sites

For model testing, the four following O'CCMoN locations were chosen: Cascade Head, Moose Mountain, Soapgrass, and Toad Creek. Each O'CCMoN location contains one forested site and one open clear-cut site (Figure 4).

**Figure 4.** EPA Oregon Crest-to-Coast Environmental Monitoring (O'CCMoN) transit sites with environmental West to East trends and details for the sites used in this study. Figure adapted from the Crest to Coast Overview document [28].

Overall, these sites span a wide range of elevations and habitat diversities between the coast and the Cascade Mountain snow zone (Table 2).

The Cascade Head open site was installed outside the Cascade Head Experimental Forest and Scenic Research Area-Forestry Sciences Laboratory (EFSRA-FSL) in the managed landscape at an elevation of 157 m (Table 2). The Cascade Head forest site is located 190 m to the northeast in a predominantly Douglas-fir forest at an elevation of 190 m (Table 2). The Moose Mountain, Soapgrass, and Toad Creek sites are positioned on the western side of the Cascades Mountain Range at increasing elevations and experience moderate to extreme weather. The Moose Mountain open site was installed within a forest clear-cut at an elevation of 668 m with the forest site located 460 m to the northeast in a predominantly Douglas-fir forest at an elevation of 658 m (Table 2). The Soapgrass open site was installed within a forest clear-cut at an elevation of 1298 m with the forest site located 1190 m to the northeast in a predominantly Douglas-fir forest also at an elevation of 1190 m (Table 2). The Toad Creek open site was installed within a forest clear-cut at an elevation of 1202 m with the forest site located 471 m to the east in a predominantly Douglas-fir forest at an elevation of 1198 m (Table 2).

The soil temperature probes were all installed in the same manner at all EPA O'CCMoN locations for both the open site and forested site. In each location, two soil temperature sensors were installed at a depth of 15 cm and 30 mm, respectively, from the mineral soil surface, i.e., just below the O-horizon [23]. The testing of the AST and AST-Solar subroutines did not involve any data collection, but rather leveraged the data collected through the EPA O'CCMoN project. These data and the details of field work can be found in the documents at the data repository [23].


**Table 2.** O'CCMoN Open and Forest Site Characteristics.

Note: All information in this table was obtained from the O'CCMoN dataset documentation, except the annual rainfall which was obtained through the PRISM 1981–2010 annual rainfall normals [29].

#### 2.5.2. AST versus AST-Solar Subroutine Setup

The AST and AST-Solar subroutines were both ran from 1 January 2005 through 31 December 2005 at a daily time step. For each site, the same O'CCMoN observed air temperature and soil moisture data were used as the data drivers for both subroutines [23]. The O'CCMoN data are measured in 30-min intervals, yet the VELMA model functions at a daily time step. To match the VELMA temporal grain, all observed O'CCMoN data were averaged to a 24-h period.

The VELMA spatial framework, per cell, contains four voxel layers; therefore, the AST and AST-Solar subroutines function under this spatial framework. Yet, the O'CCMoN dataset contains soil temperature probe data at only two depths. The AST and AST-Solar soil moisture probe depth variables for layer 1 and 2 were set to match the sensor depths of 15 cm and 30 cm [23]. Since the O'CCMoN data sites contained only two soil moisture probe depths for the sites selected, the AST and AST-Solar voxel layer three and four soil temperature results could not be evaluated and were excluded.

The AST-Solar model utilized the additional solar energy driver data. For the AST-Solar open site simulations, the CellSOLAR and MaxSOLAR variables utilized open site solar energy data. For the AST-Solar forest site simulations, the CellSOLAR variable utilized the forest site solar data, while the MaxSOLAR variable was calculated using the open site solar energy data.

The variable ReducerSOLAR is utilized as a fractional variable scaled from zero to one (Equation (12). This setup allows any solar energy units to be implemented through this method. CellSOLAR being the solar energy per location and MaxSOLAR representing the location with the most solar energy means for an open site, the CellSOLAR and MaxSOLAR values will similar if not the same. In this scenario, ReducerSOLAR will cause minimal to no reduction to the soil temperature. Conversely, forest site CellSOLAR and MaxSOLAR values will be quite different. In this scenario, ReducerSOLAR will cause a reduction in the soil temperature.

#### **3. Results**

Model results for each of the O'CCMoN sites are summarized in Table 3. Overall, the inclusion of spatially distributed solar energy improved the simulated solar temperature results. Both open and forested sites exhibit gains in accuracy, though the inclusion of spatially distributed solar energy was most beneficial for the forest sites. Below, all data are distinguished using the following attributes: site location, open versus forest environment, and soil layer. O'CCMoN sites with open versus forest

locations are listed with abbreviations under "Sites" in Table 3. Soil Layer 1 and Soil Layer 2 are referred to as SL1 and SL2, respectively.


**Table 3.** VELMA-AST and VELMA-AST3 O'CCMoN results.

The performance of the AST-Solar model at the Soapgrass site increased for soil layers 1 and 2 at both the open and forest sites compared to the AST model, but especially for soil layer 2 (Table 3). Specifically, the SGO-SL1 performance increased from a r<sup>2</sup> of 0.80 to 0.85 (Table 3; Figure 5A), while the SGO-SL2 performance increased from a r<sup>2</sup> of 0.69 to 0.90 (Table 3; Figure 5C). The SGF-SL1 performance increased from a r<sup>2</sup> of 0.69 to 0.92 (Table 3; Figure 5B), while SGF-SL2 performance increased from a r<sup>2</sup> of 0.57 to 0.89 (Table 3; Figure 5D).

The results among all sites are unique for each site, though the seasonal pattern and increased performance are similar over the year. Since the patterns are similar for each of the different sites, only the Soapgrass site results are graphically represented. Due to the 100-day PhaseLAG (spin-up) requirement of the AST model, only Julian days 101 through 365 are graphically represented.

The performance at Cascade Head (CH) increased for the AST-Solar subroutine compared with the AST subroutine for soil layers 1 and 2 for the forest site, but the performance only improved soil layer 2 of the open site. The open site soil layer 1 was the only simulation that exhibited a decrease in simulated versus observed agreement by decreasing from a r<sup>2</sup> of 0.83 to 0.77. In contrast, the CHO-SL2 performance increased from a r2 of 0.71 to 0.95. All CH14 simulations improved. The CH14-SL1 performance increased from a r<sup>2</sup> of 0.74 to 0.87, while the CH14-SL2 performance increased from a r2 of 0.71 to 0.94.

The performance of the AST-Solar subroutine compared with the AST subroutine at the Moose Mountain (MMO) increased for soil layers 1 and 2 at both the open and forest sites, particularly for soil layer 2. The MMO-SL1 performance increased from a r2 of 0.81 to 0.92 (Table 3). The MMO-SL2 performance increased from a r<sup>2</sup> of 0.67 to 0.93 (Table 3). The MMF-SL1 performance increased from a r2 of 0.89 to 0.93 (Table 3). The MMF-SL2 performance increased from a r2 of 0.70 to 0.94 (Table 3).

The performance of the AST-Solar subroutine over the AST subroutine at the Toad Creek increased for soil layers 1 and 2 at both the open and forest sites. The TCO-SL1 performance increased from a r2 of 0.82 to 0.83, while the TCO-SL2 performance increased from a r2 of 0.73 to 0.92 (Table 3). The TCF-SL1 performance increased from a r2 of 0.83 to 0.90, while the SGF-SL2 performance increased from a r2 of 0.64 to 0.89.

**Figure 5.** All panels compare observed data for the year 2005, AST model results, and AST-Solar model results. (Panel **A**) Soapgrass Open Site Soil Layer 1 temperature results. (Panel **B**) Soapgrass Forest Site Soil Layer 1 temperature results. (Panel **C**) Soapgrass Open Site Soil Layer 2 temperature results. (Panel **D**) Soapgrass Forest Site Soil Layer 2 temperature results. Simulation r2 values are calculated using Julian days 101–365. Ignoring the first 100 days prevented the AST model from being penalized for its 100-day phase lag, which only occurs during the first year of simulation.

#### **4. Discussion**

VELMA's original soil temperature model functioned well without the inclusion of spatially distributed solar energy (see Table 3; Figure 5), yet the inclusion of spatially distributed solar energy can provide significant improvements to simulated ecological processes driven by solar energy. The inclusion of local ground-level solar energy data improved VELMA's AST subroutine simulations of soil temperature from more than one perspective. First, the observed versus modeled comparisons for all sites improved, with one exception. The only exception was the CHO site, which was a landscape anomaly amongst the sites due to the station existing within a regularly maintained grass lawn at the headquarters of the research area (i.e., Cascade Head EFSRA-FSL). The TCO site received the smallest soil temperature modeling improvement with SL1 r<sup>2</sup> increasing from 0.82 to 0.83, yet the TCF site showed a significant SL2 improvement with a r2 increase of 0.73 to 0.92. The largest single layer improvement was observed at the SGF site with the SL1 r<sup>2</sup> increasing from 0.69 to 0.92 and SGF-SL2 r2 increasing from 0.57 to 0.89. It is worth reiterating that though the AST-Solar subroutine has a calibration parameter, no calibration was applied when simulating any of the sites to ensure the AST to AST-Solar estimates of soil temperature were a fair comparison of the subroutine performance.

Beyond the r<sup>2</sup> goodness-of-fit metrics, the daily variability in the modeled AST-Solar soil temperature data was reduced. This can be seen in all the graphs presented above by comparing how the AST subroutine demonstrated a repeated overestimation and underestimation of soil temperature compared to the observed data. However, the AST-Solar subroutine's noise was greatly reduced. This noise pattern in AST was due to the significant influence from the daily average air temperature driver data. Therefore, the static equations that provided oscillatory proxies for solar energy did not fully parallel the environmental phenomena of solar energy. In part, this disparity explains the resulting noisy estimations of soil temperature when solar energy was not directly included in the subroutine.

VELMA's AST soil temperature subroutine required a 100-day soil temperature simulation spin-up period. This time lag allowed for a sufficient temporal delay in the driver data utilized by Equation (6) (note that this time lag is only required for the first year of multi-year AST simulations). The purpose of the lag time was to account for the seasonal weather influence on the soil column. For the new AST-solar subroutine, this lag was no longer required due to the addition of localized daily solar energy data interacting with the existing soil moisture within layer 1.

Both the AST and AST-Solar subroutines do not model the insulation effects of snow pack [4]. Further work could include snow depth and its insulative effects on soil temperature. This would improve the soil temperature estimates in the winter, due the insulation of heat from the snow pack preventing the soil temperature from getting colder or even freezing. For the AST-Solar model, this may further improve performance with observed data (Figure 5, panels A–D) for days where snow pack persisted at Soapgrass (Figure 2) as well as goodness-of-fit metrics for the other O'CCMoN monitored sites. Similarly, further improvements in AST-Solar soil temperature estimates may be possible by including VELMA predictions of surface detritus (dead leaves and wood) and ground-level leaf biomass (proxy for leaf area index) that contribute to near-surface shading of mineral soil surfaces in open and forest sites.

The subroutine performance improvements reported here are due to the influence of local solar energy, which alters the resulting soil temperature. The prior soil temperature model predominantly utilized average air temperature as a proxy for energy, which is commonly done in watershed models. Though VELMA is a spatially distributed model, the default weather model is driven by single site location climate data. This setup resulted in homogeneous average air temperature across the simulated watershed. This is true for all forest cells and bare open prairie or forest clear-cut cells alike. By including local solar energy representation, the subsequent modeling of soil temperature was enhanced due to the improved model representation of the real world. This mainly was accomplished through the inclusion of the environmental variable solar energy that causes direct and significant influence on the phenomena soil temperature.

#### **5. Conclusions**

Watershed models are widely used to simulate the effects of land use change on the environment and the quantity and quality of hydrologic components throughout a watershed. In this paper, we demonstrated that local solar energy information improved soil temperature modeling estimates simulated by a soil temperature subroutine within a larger ecohydrological watershed model. These models were compared with observed data for soil temperatures at two depths within both open and forested environments among four observed data sites (i.e., EPA's O'CCMoN transect data) [23].

Overall, by including explicit information regarding the spatial distribution of solar energy across a landscape, watershed models can better capture the spatiotemporal variations of soil temperature in both forested and open sites. Therefore, researchers that utilize spatially distributed or semi-distributed mechanistic watershed models should consider incorporating spatially explicit solar energy models (e.g., Penumbra [30,31]) or other spatially heterogeneous descriptions of ground-level solar energy to better represent energy exchange at the surface. This is especially true when modeling discrete landcover types such as forested, open, water, and agricultural cover and when modeling the impacts of riparian shading on soil temperature and stream temperature, as well as the effect of solar energy on fish habitat.

Finally, while we presented improvements of a soil temperature subroutine within an ecohydrological model, other subroutines can also benefit by the inclusion of spatiotemporal representations of ground-level solar energy. The integration of local solar energy information with watershed models and all their subroutines could potentially benefit several processes, such as snow melt, water temperature, and plant growth via photosynthesis. Integration of spatially explicit ground-level solar energy models with environmental models can provide dynamic feedbacks between other environmental processes, such as tree growth and shade. As tree growth is simulated within watershed models, their heights could be transferred back to the ground-level solar energy model in a dynamic mechanism, which then would alter the amount of solar energy that is intercepted by the canopy and does not reach the ground. This dynamic integrated modeling approach could be extremely beneficial for looking at the long-term effects of planting riparian buffers and determining the duration required for stream temperatures to be reduced by some threshold. Because solar energy is amongst the most impactful environmental variables in natural and managed ecosystems, further investigations would greatly benefit from the application of watershed-scale models that are dynamically coupled with spatially-explicit solar energy models.

**Author Contributions:** J.J.H., B.L.B., R.B.M., R.E.K. and J.J.G. designed the project and modeling approach; J.J.H. and B.L.B. performed, and R.B.M., R.E.K., A.F.B., K.S.D., and P.P.P. assisted with, the model set up and simulations; R.S.W performed all field data curation; J.J.H., B.L.B., R.E.K., R.B.M., J.J.G. and R.S.W. analyzed the results and developed the manuscript's structure; all authors contributed to writing the paper.

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

**Conflicts of Interest:** The authors declare no conflict of interest. Findings and conclusions in this article are those of the authors and do not necessarily reflect or represent the views or policies of the U.S. Environmental Protection Agency. Any mention of trade names, products, or services does not imply an endorsement by the U.S. Government or the U.S. Environmental Protection Agency. The EPA does not endorse any commercial products, services, or enterprises.

#### **References**


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

#### *Article*

## **Issues of Meander Development: Land Degradation or Ecological Value? The Example of the Sajó River, Hungary**

#### **László Bertalan 1,\*, Tibor József Novák 2, Zoltán Németh 3, Jesús Rodrigo-Comino 4,5, Ádám Kertész <sup>6</sup> and Szilárd Szabó 1,6**


Received: 29 September 2018; Accepted: 6 November 2018; Published: 9 November 2018

**Abstract:** The extensive destruction of arable lands by the process of lateral bank erosion is a major issue for the alluvial meandering type of rivers all around the world. Nowadays, land managers, stakeholders, and scientists are discussing how this process affects the surrounding landscapes. Usually, due to a land mismanagement of agroforestry activities or urbanization plans, river regulations are designed to reduce anthropogenic impacts such as bank erosion, but many of these regulations resulted in a degradation of habitat diversity. Regardless, there is a lack of information about the possible positive effects of meandering from the ecological point of view. Therefore, the main aim of this study was to investigate a 2.12 km long meandering sub-reach of Sajó River, Hungary, in order to evaluate whether the process of meander development can be evaluated as a land degradation processes or whether it can enhance ecological conservation and sustainability. To achieve this goal, an archive of aerial imagery and UAV (Unmanned Aerial Vehicle)-surveys was used to provide a consistent database for a landscape metrics-based analysis to reveal changes in landscape ecological dynamics. Moreover, an ornithological survey was also carried out to assess the composition and diversity of the avifauna. The forest cover was developed in a remarkable pattern, finding a linear relationship between its rate and channel sinuosity. An increase in forest areas did not enhance the rate of landscape diversity since only its distribution became more compact. Eroding riverbanks provided important nesting sites for colonies of protected and regionally declining migratory bird species such as the sand martin. We revealed that almost 70 years were enough to gain a new habitat system along the river as the linear channel formed to a meandering and more natural state.

**Keywords:** bank erosion; landscape metrics; diversity; Sajó River; UAV

#### **1. Introduction**

Alluvial rivers represent dynamic landforms of the watersheds all over the world [1]. Under minimal regulation, these rivers can generate diverse landscapes by shifting back and forth along

their floodplain [2,3]. In meandering rivers, the channel flow undercuts the outside banks resulting in seepage outflow or mass failure processes i.e., slab failures, rotational slides and sloughings or slump blocks eroding into the water body [4,5]. These materials can form point bar accretion at the inside banks downstream due to the lower current velocity at the inner side, as opposed to the outer side of the bends [6]. The cross-sections of meandering riverbeds are complex showing asymmetric bathymetric and flow patterns, which is reflected in habitat diversity as well [7]. Point bars are specific fluvial geomorphic features inside a meandering river bend [8]. They are deposited from finer sediments by the low energy parts of the river along the inside bank driven by recirculation zones [9]. The ridge and swale topography produce undulating bar surfaces [10] and the long-inundated periods of the depressions allow pioneer bush and tree species colonizing the fine-grained surfaces. Intensive channel migration maintains ideal conditions for vegetation succession processes [11], since the bend migrates forward through lateral aggradation, the colonization process increases [12]. Generally, the maxima of lateral bank erosion rates are concentrated downstream of the point bar symmetry-axis [11] resulting elongated, skewed or compound meander bends [13].

River channels and their surrounding floodplains enhance landscape evolution and the diversification of environments [14]. Quasi-natural rivers without extensive channelization, bank protection and embankments, are susceptible to rapid changes in their hydro-geomorphological features. One of the main reasons for this is that vegetation colonizing riparian zones are affected by frequent flood-disturbances or even severe droughts, therefore, they have adapted to variable water levels. These plants have a major influence on the initiation of geomorphic and hydraulic processes of the channel and floodplain as they regularly stabilize organic matter from sediment fluxes [15]. The majority of these species need sunlight for growth that is not available under dense canopy levels; therefore, they colonize freshly formed open spaces, and the resulting patchy vegetation mosaics well-represent the habitat diversity established through flooding, scouring and sedimentation [16]. Furthermore, the resistance properties of the vegetation on channel flow is recognized as a potential factor in the mitigation of land degradation processes [17,18]. River corridors, as linear features of fluvial landscapes, are such integrated ecological systems that connect identical landscape elements. Due to the dynamic interactions within climatic factors, catchment geology, relief, inundation and nutrients, the ecological turnover can be outstanding in these features [12]. The aquatic and terrestrial vegetation patches provide heterogeneous and diverse habitat for fish, water birds and macroinvertebrates [19]. The initial colonizers maintain their stands, preventing further recruitment of vegetation beneath their canopy level [20,21]. Furthermore, vegetation patches substantially block and divert river flow around their canopy, and flow velocities decrease along the vegetation patch and increase in the surrounding channel area; thus, they determine the future occupancy patterns of trees along the point bars [22]. The succession of even ephemerally emerging bar surfaces provides suitable foraging and nesting sites for water birds. Moreover, due to their isolated nature, these sites are often free from human disturbance as well [19,23].

Lateral bank erosion of meandering rivers is responsible for extensive destruction of arable lands and usually threatens human environments [24,25]. These phenomena, accompanied by riparian deforestation and serious flood inundation, provide the necessity for the protection of most European river networks [26–28]. Regulated river channels have less cross-sectional diversity from both geomorphological and ecological point of view. In the past years, several studies revealed that channelization and bank protection works are responsible for an extensive, global-scale habitat and river ecosystem degradation [29–32]. It is also proven that river regulation impacted fish and macroinvertebrates negatively due to the degradation of habitat heterogeneity [33]. Recent studies found 're-meandering' as an effective practice for the restoration of habitat diversification at shortand medium-terms [34], since lateral erosion and the channel migration maintain sediment supply for vegetation colonization [35]. However, long-term ecological consequences of re-meandering have not been well-studied yet [36]. Only a few studies found positive trends in diversification but many

restoration projects were not able to demonstrate significant differences in biodiversity between the before and after periods of the restoration project over short (5 to 10 years) time periods [9,37,38].

In Hungary, the majority of rivers had been regulated and channelized, starting in the 19th century [39–42]. However, few medium-size rivers of the Tisza River Basin such as the Sajó River remained in a quasi-natural state since economic issues had prevented the total channelization [43]. Along with several reaches of Sajó River, extensive lateral erosion threatens the agricultural areas. However, there is a lack of information about this vital issue, which could be included as key information for land management plans. Thus, the main aims of this research were (i) to investigate the geomorphological development and effects of bank erosion along a meandering sub-reach of Sajó River; and, (ii) to assess the possible impacts on ecological diversity in case there will not be further interventions on channel morphology. To achieve these goals, we performed a GIS (Geographic Information System)-analysis of three consecutive meandering bends over 10 periods between 1952 and 2017 based on archive aerial imagery and UAV-surveys. Moreover, an ornithological survey was also carried out to assess the composition and diversity of the avifauna.

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

#### *2.1. Study Area*

Sajó River (or Slaná in Slovakia) is a transboundary river of Slovakia and Hungary having a total length of 229 km (124 km in Hungary). It is the main tributary of Tisza River before it reaches the Great Hungarian Plain. The river catchment is situated at the Eastern-Carpathians with a total area of 5545 km2. The stream gradient is much higher in Slovakian territory, then the river becomes alluvial downstream from the Hungarian state border. The mean discharge of the river is around 24 m3/s. The total average of suspended sediment load varies between 828,000 m3/y and 1,927,000 m3/y [44]. The Hungarian reach of Sajó River is mainly the alluvial meandering-type with a total sinuosity of 1.78.

Although extensive river regulation plans had been established in the early 20th century to facilitate shipping towards the Carpathians, eventually, most of the works had to be cancelled due to the economic issues associated with World Wars [45]. Only a few minor river management works had been carried out mainly around industrial areas. Most of these are artificial cutoffs, small groynes, and another bank protection. Even though their spatial distribution is broad (58.3%) [43,46] along the Hungarian reach of Sajó River, this is one of the least regulated rivers in the country. The abovementioned geographical settings and the low rate of human intervention lead to the situation of lateral bank erosion that can reach the 4–7 m/y rate in several sub-reaches [46,47].

The Hungarian National Ecological Network (Phase 2) had been established between 1999–2001 according to the legislation of the Pan-European Ecological Network (PEEN). The components of the network provide the following landscape element categories: well-known core-areas, ecological corridors, buffer zones and restoration areas [48]. Core areas provide the main habitats and genetic reserves while the strip-like ecological corridors serve as continuous habitats or chains linking the smaller and larger habitat patches together [49]. This study focuses on a selected 2.12 km long sub-reach of the Sajó River consisting of three consecutive river bends near the town of Nagycsécs, Hungary (Figure 1). The further calculations were performed on this 68.4 ha large rectangular area around the river channel. This sub-reach of Sajó River can be considered as a free-forming meandering type located between major river engineering works both upstream and downstream. The study area as part of the Sajó River floodplain belongs to the PEEN category of ecological corridors between two main Natura2000 areas.

**Figure 1.** Overview of the study area. (**A**) Location of the study area in Europe; (**B**) Location of the main rivers in the region of study area; (**C**) Detailed overview of the selected sub-reach of Sajó River.

We observed that the land mosaic of the study area was represented by a patch-complex consisting of bare surfaces, perennial grasslands, forests and bushes, the river channel, arable lands, and settlements. Bare plots, grasslands, and arboreal vegetation were considered as patches connected to each other by their own successive development controlled by the river channel changes and flood dynamic. On the other, arable lands and settlements were considered as mainly human-controlled elements of this structure. Bare surfaces meant patches not completely in a lack of vegetation, but with a low cover. Plant species of these sites are mainly disturbance-tolerant species such as *Chenopodium album*, *Chenopodium ambrosioides*, *Bidens tripartitus*, at their wet riverside edges *Polygonum hydropier*, *Polygonum minus*, *Polygonum mite*, *Rorippa amphibia*, *Rumex crispus*, and *Rumex obtusifolius*. Perennial grasslands are dominated by diverse grass species like *Phalaris arundinacea*, *Agrostis stolonifera*, *Alopecurus pratensis*, and *Agropyron repens*, varying by land use and duration of floods. In forest and bush patches grow *Salix purpurea*, *Salix triandra*, *Salix alba*, and *Populus* spp.

#### *2.2. Datasets*

This study aims to develop a spatiotemporal analysis on the river channel development and the ecological diversity based on a set of aerial imagery in 10 different time periods (Table 1). A set of black/white military-based historical aerial photographs of the study were available from the year of 1952, 1956, 1975 and 1988 given by the Hungarian Military History Museum. The archive aerial imagery was scanned in 600 dpi resolution and then orthorectified using the ERDAS Imagine software. However, orthophotographs of 2011 were used as a reference (datum: Hungarian HD72/EOV), but we included

Topographical maps of 1980 as well since there were few field objects on the aerial imagery between 1952 and 1988, that w not possible to identify on the orthophotographs of 2011. During the process, 15–20 ground control points (GCP) were used for each image to provide accurate georeferencing. The RMSE (Root Mean Square Error) values varied between 3.4 to 6.7 m with an average of 4.8 m. We purchased digital colour orthophotographs of 2000, 2005 and 2011. The scale range of the aerial imagery and orthophotographs were found between 1:7000 and 1:12,000. After 2011, there were no official orthophotographs available from the study area. Nowadays, UAVs offer a valuable solution to produce high-resolution aerial imagery from a few km2 large areas [50,51]; thus they are widely used for mapping wetland areas, especially for disaster management [52,53]. UAV-based surveys were conducted in 2015 by using DJI Phantom drones in order to provide orthophotographs for 2015, 2016 and 2017. Each flight was performed at low-flow conditions while, at least, 20 GCPs were measured by a Stonex RTK-GPS system. UAV image acquisition was performed at 150 m A.G.L. with 75% frontlap and sidelap between flight paths to provide a ground resolution of 0.07 to 0.09 m. Agisoft Photoscan software was used for the photogrammetric processing and the creation of orthophotographs with an overall accuracy of around 0.05 m.


**Table 1.** Basic parameters of aerial imagery and orthophotographs used in this study.

#### *2.3. Indices of River Channel Development*

In order to quantify the extent of degradation by the lateral bank erosion, overlays of pairs (Figure 2) related to consecutive time periods of the river channel polygons were composed [3]. In this research, we consider that understanding the area of accretion also plays a key role in the case of ecosystem diversity analysis. These two variables were derived at the same time. The crosshatched area represents that part of the floodplain where the two-channel planforms overlap and it appears that the channel did not change position, without any pronounced erosion or accretion.

**Figure 2.** Methodology for calculating areas of erosion and accretion.

The value of sinuosity index (SI) was calculated as the ratio of channel length to valley length [54]; consequently, not just the areal but morphological trend of river channel development could be

assessed too. We determined the channel length as the length of channel centerline while valley length was measured as the distance between the point where Sajó River crosses the Hungarian border and the point of the Tisza River confluence. We calculated total SI values for each investigated time-period. We determined the mean channel migration rate based on polygons drawn by the overlapping centerlines in each consecutive time periods. The mean lateral channel shift is a ratio of polygon area and the half of the polygon perimeter [55]. Regarding erosion/accretion and channel shift, we normalized the values by the number of years covering each different periods.

We identified the meander parameters e.g., chord length (straight line distance between two inflexion points; but this value is not equal with meander wavelength, that is the straight line distance of two consecutive meander apexes of the same side of the river banks), amplitude, the width-normalized radius of curvature (R/w). Channel widths were also determined in every 20 m along the centerlines of each time periods; therefore, a detailed mean channel width was given at the end.

#### *2.4. Ornithological Survey*

The avifauna of the study area via transect and point count surveys during the breeding season of 2018 was assessed. We selected two 200 m long line transects traversing the recently deposited areas as well as an observation point along an eroding section of the river to be able to detect birds along the full length of the study area. The transects were 500 m apart, while the point was 225 m and 580 m from the transects. We conducted the survey on 26 June between 10:00 a.m. and 12:40 p.m. Each transect survey lasted for 30 min and the point count for 20 min during which the observer recorded all species seen or heard in the vicinity of the river channel and the surrounding floodplain that are potentially using the habitat (i.e., not just flying over it at a high altitude).

We also quantified the size of the sand martin colonies of the eroding river banks of the study area from photo-mosaics. We applied close-range terrestrial photogrammetry to compose the mosaic view of the eroding river bank. For this purpose, a Nikon D5300 camera with a 70–300 mm tele-objective lens was used to capture the necessary photo-sequences of 218 images from tripod stands from the opposite side of the river. The datasets were processed in Agisoft Photoscan software. We printed coded targets of the software in A4 size in order to place them on the bank edges. We measured their precise coordinates by Stonex RTK-GPS system in order to provide a reliable scaling of the mosaic in Agisoft Photoscan.

#### *2.5. Landscape Metrics*

Class and landscape level landscape metrics to reveal the ecological features and processes along the changing river channel were applied. Initially, landscape patches were manually vectorized with a scale of 1:1000 in ArcGIS 10.3 on the aerial imagery in each investigated time-period. During the process, we have delineated the following land use categories: forests and bushes, grasslands, arable lands, bare point bar surfaces, settlement and the river channel itself (Figure 3). Aerial photographs of 1952, 1956, 1975 and 1988 are black and white ones and sometimes the quality is far from good, however, these photographs ensure a larger time range to monitor the changes. Land cover classes had to be consistent; therefore, we merged the bushes and trees of forests. We decided to apply this approach considering the four stage successional model of Corenblit [56]. In our case bare surfaces represents the vegetation recruitment phase, grasslands are in phase of vegetation establishment and succession continues with development of bushes and forests, which consists the same spectrum of species but in different age and development phase. The river channel itself is the main factor for diaspore dispersion. Arable lands and settlements are totally controlled by the society, until the migration of the river channel does not alter this situation and shift them to a previous phase of succession.

**Figure 3.** Aerial overview of the study area with the different land use categories (1—forests and bushes; 2—grasslands; 3—arable lands; 4—bar surfaces; 5—river channel). (**A**) Oblique aerial photo by L.B., 11 June 2017; (**B**) Black & White archive aerial photo from 1975; (**C**) UAV-based orthophoto from 2017.

In order to quantify and evaluate the direct changes in land cover and the transformation of land patches related to vegetation succession, channel migration or human agricultural management, we calculated a confusion matrix based on the vectorized land cover files in Idrisi software. In these tables, the diagonal values related to a pair of a given land cover category represents the proportion of area where no changes were found. The columns represents the initial land cover category while the rows will define the land cover categories that the given column category were transformed into. The values represents their proportion from the total area.

We concentrated on the metrics which reflect the diversity of the floodplain and its environments:

• Patch Density (PD): we calculated the index on landscape level as the number of patches per unit area (Equation (1)).

$$PD = \frac{N}{A}(10,000)(100)\tag{1}$$

where *N* is the number of patches in the landscape, and *A* is the total area (m2), and the outcome is expressed in number. per 100 ha−<sup>1</sup> [57].

• Interspersion and Juxtaposition Index (IJI): we calculated the index on landscape level as the function of observed interspersion and the maximum possible interspersion for the given number of patch types (Equation (2)).

$$\text{LII} = \frac{-\sum\_{i=1}^{m} \sum\_{k=i+1}^{m} [(E\_{ik}) \times \ln(E\_{ik})]}{\ln\left(\frac{m(m-1)}{2}\right)} \tag{2}$$

where *Eik*: the total edge between patch types *i* and *k*; *m*: number of patch types. IJI is expressed in per cent (0–100); low values indicate unevenly distributed or isolated patches in the area, the largest value is acquired when all patch types have a common edge with all possible other patch types [57,58].

• Shannon's Diversity Index (SHDI): We calculated the index on landscape level as the sum of the proportional abundance of each patch type multiplied by that proportions (Equation (3)).

$$\text{SHDI} = -\sum\_{i=1}^{m} (P\_i.InP\_i) \tag{3}$$

where *Pi* is the proportion of the landscape occupied by the patch type *i*. When the landscape is occupied by only one patch SHDI = 0 and increases with the number of patch types and their proportional area without upper limit [57,59].

• Class Area of the forests (CA\_F): we considered the forests' area as the indicator of landscape change (in the initial phase, in 1956, there were only a few small patches and later, with the river bed development, the area relevantly increased). The index was calculated as the proportion of the forest land cover type and the total area and was expressed in per cent.

#### *2.6. Statistical Analysis*

After the calculation of descriptive statistics, the connection between the indices of river channel development (Sinuosity—SI, Area of erosion—Er, Area of Accretion—Ar) and the landscape metrics was also investigated. Regression analysis and Principal Component Analysis (PCA) were applied. While the regression analysis pointed on the changes directly with the given variables and provided valuable information on the temporal features and breakpoints of the surface development, PCA helped to identify the years when both the landscape and the river bed had deterministic relations.

Standardized PCA was conducted using the correlation matrix with Varimax rotation to gain orthogonal axes, i.e., non-correlating principal components (PCs) [60,61]. PD, CA\_F, IJI, and SHDI landscape metrics and the Er, Acc and SI indices of riverbed development were involved. Thus, it became possible to identify both the cross-connections among the variables and to visualize the dates of different stages of surface development on a biplot diagram. Model fit was controlled with the Root Mean Square Residual (RMSR), which is determined using the residuals of the original correlation matrix and the estimation of the PCA [62]. RMSR values of <0.05 indicates very good fit [63].

We applied the Jonckheere-Terpra test to reveal whether there was a trend in river channel change over the studied period. Statistical analysis was performed using the software Past 3.19 [64] and R 3.5 [65] by applying the psych [66] and GPArotation [67] packages; furthermore, the lattice [68], the clinfun [69] and ggplot2 [70] packages were applied for the data visualization.

#### **3. Results**

#### *3.1. Changes in Land Cover and Channel Morphology*

Detailed land cover changes based on the vectorization of the available aerial imagery are shown in Figure 4. After a visual interpretation of the figure, it is clearly visible that significant changes occurred over the investigated periods (65 years). The channel developed at a remarkable rate developing three main meandering bends with high sinuosity.

According to the results of the normalized erosion and accretion rates, the first period between 1952 and 1956 showed the second highest bank erosion (0.85 ha·year<sup>−</sup>1) rate (Figure 5A). Bank erosion activity radically decreased with three times by 1975. It was followed by a gentle increase between 1975–1980 but after that, it decreased again and reached the lowest bank erosion rate (0.16 ha·year<sup>−</sup>1) in 2005. A notable increase started in 2010 with an outstanding maximum value of 1.12 ha·year−<sup>1</sup> in the period of 2015–2016. Except for the first period (1952–1956), accretion rates followed closely the erosion rates; moreover, from the period of 2000–2005, except the above-mentioned extremely erosive year between 2015 and 2016, they exceeded their contribution over the eroded areas in the channel development. Lateral migration rates (Figure 5B) primarily followed the trend of erosion rates as it started from a relatively high rate (5.3 m) of channel shift and decreased by 1975. Its minimum

(0.8 m) was also found in 2005 then similarly started to increase, however, its peak of 2016 was not as outstanding, due to its erosion rates.

**Figure 5.** (**A**) Time series of normalized accretion and erosion rates in the studied period; (**B**) Time series of mean lateral channel shift rates in the studied period.

The parameters of meander evolution are summarized in Figure 6. The chord length (Figure 6A) increased in Bend 1 from 460 m to 634 m until 1988 but it was followed by a slight decrease. However, Bend 2 showed the opposite trend as it decreased to 299 m by 1988, but then the chord length increased intensely with almost 200 m by 2017. Bend 3 decreased monotonously and the total change was almost 350 m from 1952 to 2017. Amplitude (Figure 6B) increased in all the bends but the expansion of Bend 1

69

was almost two times higher (+127.49 m) than the others. The change in width-normalized radius of curvature (Figure 6C) showed a similar increase as it was in chord length between 1952 and 1988 but it was also found at Bend 3 as well. It was followed by an intensive decrease, especially in Bend 3.

**Figure 6.** Parameters of meander evolution and channel shift (**A**) chord length (m); (**B**) bend amplitude (m); (**C**) width-normalized radius of curvature (m).

We justified a significant decreasing trend in case of river channel width (Figure 7) as a function of time (J-T statistic = 80298; *p* = 0.0004) and there was a negative correlation between the mean channel width and sinuosity (*r* = −0.93, *p* < 0.001).

**Figure 7.** Changes in the mean channel width.

Land cover changes can be divided into two groups: (1) settlement, river channel and point bars where the changes were minimal (the area was more or less constant but the patches had spatial changes); (2) arable land, grassland and forests where the changes were relevant (these patch types turned into another, at least, partly). The area of grasslands decreased transforming into arable land and forest (Figure 8). The extent of the forests and bushes increased by 25% since the share of grasslands showed a decrease in similar rate and bare surfaces were not changed considerably. The extent of human-controlled arable lands was increased at the beginning but slightly decreased in the last decades. Extents of the river channel and settlement were not changed considerably.

**Figure 8.** Changes in land cover over the studied period by patch types.

The detailed contingency tables of land cover changes between each consecutive time periods were summarized and available in Supplementary File S1. In Table 2, we presented the changes of selected pairs of land use category transformations that showed higher conversion proportions. The first four rows represent the successional phases that follow the channel migration: The second highest transformation rates were connected to the transformation from former river channel to bar surfaces. It was followed by colonization of the bar surfaces by grasslands where the transformation rates were found also lower and totally stopped by 2005. However, generally, direct transformation from bar surfaces to forests showed some higher values but it could have been expressed by the longer time periods covered. The highest transformation rates were found in the transformation from grasslands to forests. The last four rows represents transformation categories that were controlled by the channel migration and bank retreat except the changes from grasslands to arable lands. This type also resulted in an outstanding proportion (0.326) between 1975 and 1988 and it radically decreased afterwards.



AL—Arable lands, BS—Bar surfaces, F—Forests, G—Grasslands, RC—River channel.

A linear relationship between the increase of forest areas and the SI (*R*<sup>2</sup> = 0.94, *p* < 0.001; Figure 9A) was revealed. At the first date of the study period (1956), only a few patches of the forest were observed, which reached a 25% proportion till 2017; besides, SI indicated a straight river channel section in 1956 and it became a complex form up to 2017. Thus, the forest areas increased directly as the river channel transformed and appropriate areas developed in the floodplain.

**Figure 9.** Sinuosity as independent factor of land cover changes and landscape diversity in the study area. (**A**) The relationship between the SI (Sinuosity Index) and CAF (Class area of Forests); (**B**) Relationship between the SI (Sinuosity Index) and SHDI (Shannon's Diversity Index) metrics.

There was a significant connection between the SI and SHDI variables (*R*<sup>2</sup> = 0.93, *p* < 0.001). However, unlike in the case of forest areas, the connection was not linear (Figure 7B); a second-order polynomial regression indicated that there was a change in the trend after 2011: next years (2015–2017) were outside of the linear trend and the curve showed saturation; moreover, there was a small decrease in the diversity in 2017.

PCA explained the 94% of the total variance. Three principal components (PC) were justified by the RMRS (0.03). PC1 correlated with the Er (Area of Erosion), Acc (Area of Accumulation) and PD (Patch Density) explaining 40% of the total variance. PC2 correlated with the SI, SHDI, and CAF and accounted for 37% of the total variance, while PC3 was formed by solely the IJI having 18% in the explanation of the total variance. Acc and Er correlated with the PD and SI with the SHDI and CAF.

The ordination diagram showed that years formed two distinct groups differentiated by the Er, Acc and PD indices (Figure 10). The first group was situated in the negative part of the diagram, and the second, with a larger range, in the positive part. Considering the vertical axis (PC2), the distribution followed a monotonous increasing trend between 1956 and 1988 but it became sparse after 2000 with 2005 having the lowest value and after another increase, it started to decrease from 2015.

#### *3.2. Avifauna*

Since the biodiversity can be affected by the different land cover and morphological river changes, we surveyed the avifauna as a possible bioindicator. Result of the snapshot faunistic survey is shown in Table 3. A total of 26 species were detected breeding or using the habitat during the breeding season, of which 23 are protected under the Hungarian law including the strictly protected European bee-eater (*Merops apiaster*) [71].

**Figure 10.** Biplot of the PCA model conducted with the landscape metrics and river channel development indices (dashed arrows: involved variables).


**Table 3.** List of bird species seen and/or heard during the censuses listed in taxonomic order.

The nest cavities of both the bee-eater and the steeply declining Sand martin (*Riparia riparia*) was found along the eroding section of the Sajó River bank. Approximately, 80–100 pairs of sand martins were feeding young at the time of the survey at the focal section of the river. However, based on the number of visible nest cavities (Figure 11), the two colonies could have consisted of a total of 446 breeding pairs.

**Figure 11.** Nest cavities along the eroding river banks of Sajó River at the study area. (**A**) The overview of the cavities; (**B**) Examples of cavities along the recently eroded slump blocks; (**C**) Red circles indicates the cavities identified on one part of the mosaic.

#### **4. Discussion**

Fluvial morphological changes are able to generate intensive modifications, which can affect hydrological processes and, even, biodiversity. One clear example is the Sajó River, which was assessed in this research. Sajó River itself is regulated only in particular reaches, but its morphological evolution is affected by the Tisza River. Sajó River as a tributary joins Tisza River and the study area is located only about 20 km upstream from the confluence, which is strongly affects the sediment transport and channel development of the studied Sajó River sub-reach. The Tisza was regulated during the second half of the 19th century, and as a result of these regulations, its riverbed was incised for its straightened channel and increased energy [42,43]. This process could be also responsible for the morphological changes of Sajó River. The characteristics of Tisza channel development basically changed after the construction dams on Tisza River, both on upstream (1954, Tiszalök) and on downstream (1978, Kisköre) sections around the firth of Sajó. The establishment of these dams could cause a reduction of sediment supply [54] which generally leads to further channel incision. However, we determined the narrowing (Figure 7), and therefore, similarly to Tisza River, a possible incision, of the Sajó River but its channel started to develop in a different way, namely increased meandering started instead of a former incision, which has diverse effects on land cover, conservational value and productivity on surrounding landscape [72].

Morphological changes allow registering high rates of bank erosion, which in similar investigations in other European rivers, was also calculated [3,73,74]; however, the extent of arable lands lost by lateral erosion is not an outstanding value compared to other rivers with similar geomorphological patterns [75]. In this way, it would be a great opportunity to include in future investigations the correlation between these morphological alterations, the changes in the biodiversity and the bank erosion rates.

It is notable that the accretion of the river bank mainly followed the rate of eroding outer banks; therefore, new bare point bar surfaces had been developed in a short time period. Although these surfaces are not valuable ecologically, at the beginning phase, they represent potential future habitats [74]; our study pointed out that in case the vegetation can rapidly occupy the area, then the process can increase the extent of habitat in the ecological corridors. Moreover, the studied sub-reach is situated between two Natura2000 areas, thus, channel migration can enhance the habitat diversity and species connectivity between these sub-groups of the ecological network, that strengthen the role of an ecological corridor, as a possible compensation of the land degradation processes.

The majority of previous studies were focusing on the extensive negative effects of river bank erosion i.e., spreading pollutants from upstream reservoirs or remobilizing heavy metal contaminants [76–78] but only a few discuss their ecological consequences. In some cases, under different climate conditions, the lateral shifts of river channels rapidly decrease the biodiversity of the neighbouring flora and fauna [79] but in this case, it affects only agricultural areas. In our study, we demonstrated an opposite dynamic: How the process can be valuable. Regarding the changes in landscape metrics, as SI values increased, the forest areas (CA-F) also increased. SI reflected how the planform of the river channel was developed, and the positive connection with the forest areas indicate that instead of the arable lands, a natural afforestation process was initiated. Small grassland patches can be merged into forest cover without appropriate management [80] but in this case, the increase of SI provided the criteria of the forest cover extension. The newly developed point bars were occupied by plants (firstly with herbaceous plants) and later became forested following a successional process. Diversity did not follow the linear trend because the area of the forests reached a threshold when the further increase did not increase the diversity of the mosaic of land cover patches (Figure 9). This is the explanation of the two kinds of the relationship of the SI with the landscape metrics. The analysis on the conversion matrix between land cover categories showed higher values on the changes related to vegetation succession phases while the lower values were found in connection with the channel migration. However, we found outstanding transformation from grasslands to arable lands between 1975 and 1980 but this can be clearly explained by human intervention on the land management.

Multivariate analysis revealed that 1956–1988 period had rather similar features considering the erosion-accretion processes and the PD, and there was another group of years between 2000 and 2017 that can be distinguished according to the remaining variables. However, considering the results of PCA, SI, SHDI, and CA\_F variables differentiated the points a sparse distribution: As the erosion and accretion processes increased or decreased, or were larger or smaller compared to each other, the forest areas and the diversity also reacted differently. This process became the major regulator of the ecological process from 1988 to 2000 when the river channel development reached SI = 1.4 and in 2015 when it reached SI = 1.6, the pace of the lateral movement slowed down. At the beginning of meander development, the chord length remains stable and an extension, with increasing amplitude, occurs [11]. We observed similar changes along the studied Sajó River reach since the chord lengths showed only minor changes until 1975 from where both increase (bend 2) and decrease (bend 1, 3) took place. According to [13] the decreasing normalized curvature (R/w) accelerates the channel migration as it was also found (Figures 5b and 6c) in the studied Sajó River reach. The proportion of the forests increased from 1 to 25%, and, what is important, these areas were constant (if a patch type transformed into forest, most of the area remained forest in the consecutive years, too), while the grasslands, arable lands, and especially the point bars changed in spatial terms (changes can occur back and forth, or the formation and diminishing is a common process).

The development of successional riparian and fluvial forests on recently deposited point bars significantly increases biodiversity locally [81] but more importantly, it facilitates the movement of organisms across an otherwise cultivated and homogenous landscape, thus, allows the maintenance of gene flow between meta-populations. Indeed, the classification of the study area is "ecological corridor" within the National Ecological Network, which connects ecologically important core areas along the river. The avifauna of the study area is similar to communities detected along with other reaches of the river over two decades ago [82] suggesting that the newly formed habitats are quickly colonized by protected species from nearby areas.

The eroding banks provide important nesting sites for colonies of protected and regionally declining migratory bird species such as the sand martin and the European bee-eater [83] further increasing the ecological importance of the new habitats created by this dynamic river. In fact, an effective way to improve the availability of nesting sites and facilitate the recolonization of an area by sand martins is the removal of bank erosion control projects from heavily regulated river channels [84]. This may be particularly important for maintaining viable populations of sand martins in agricultural landscapes as these birds usually do not reuse their nests from previous years to avoid the costs of heavy infestation with parasites [85,86], rather select nesting sites along river banks and sand quarries where periodically renewing vertical surfaces are available at the beginning of the nesting season. Eroding river banks at our study area provide just that, the continuous renewal of this critical resource, the nesting site, which often limits the distribution of both sand martins and bee-eaters [71,85].

The increase in the extent of the most stable successive landscape element (forests) and a decrease of elements in middle-phase of succession (grasslands) suggest that the land-mosaic alteration during the investigated time resulted in higher diversity and stability. This process is important regarding laterally active channels; therefore, identification of the relationship between morphological changes in river channel and landscape evolution is vital [2,56]. Changes of arable lands could not be used as an indication for river development since it is controlled by human actions related to agricultural management. As a main controlling factor for successive landscape elements, the altered dynamic of meandering could be evaluated. Even if no relevant human interventions were present within the study area, the changes in meandering-dynamic can be highly influenced by the anthropogenic changes in the flood-dynamic of the main stem Tisza River. Decreasing agricultural productivity due to bank failures and channel shifting seems to be balanced by increasing habitat diversity by recent point bar deposits, succession areas, and bushes. Considering the increasing share of cultivated areas during last decades within a floodplain, which is part of a national ecological network, positive effects of the meandering of the river exceeds the negative effects of land loss and lateral erosion.

Another important topic to be assessed in the future would be the connectivity processes [87]. It is clear that the studied river and surrounding landscapes are connected by different processes such as nutrient transport or sediment mobilization, both of these are also conditioned by the vegetation colonization [88–90]. The application of modelling techniques and connectivity indexes allow us detecting how meandering is influenced by other dynamic fluxes such as agricultural fields [91] or urban areas [92] and at which level. Therefore, undoubtedly, understanding connectivity processes in the Sajó river would help land plan managers and stakeholders to design correct and sustainable soil erosion control measures and water conservation practices, as other authors also confirmed in other degraded areas [93,94].

#### **5. Conclusions**

River channel development is a natural process, which usually considered to have a negative effect due to the damages caused to agriculture or infrastructure. Our hypothesis was that, besides the detrimental phase, there is also a positive effect because of the newly developed habitats on the opposite side of the river. We revealed that 65 years were enough to gain a new habitat system along the river as the linear channel formed into a meandering and more natural state. At the beginning phase, there were only a few patches of forests and the matrix had been dominated by the surrounding arable lands, while nowadays the forests have a major role in the landscape mosaic. There was a linear relationship between the sinuosity and the class area of the forests; i.e., the more developed the meanders were, the more forest patches appeared in the area. However, Shannon's diversity did not follow a linear trend with sinuosity and instead of reaching its maximum it showed a polynomial trend; i.e., the areal increase of forests did not increase the landscape diversity as the patches became dominant, compact and connected. Consequently, this was an advantageous process from ecological aspects. Although the plant species were not of significant conservation value (mostly pioneers and weeds) but provided habitats for several protected bird species. Besides, the eroding side of the riverbed also serves as a nesting place for birds, too. Accordingly, we emphasize the positive effects of the erosion and accretion processes, as nature conservation benefits from the new geomorphological forms of river channel development.

**Supplementary Materials:** The Supplementary Materials are available online at http://www.mdpi.com/2073- 4441/10/11/1613/s1.

**Author Contributions:** L.B. designed and carried out this research; S.S. and Á.K. supervised and instructed this research; L.B., S.S., T.J.N., and Z.N. prepared and analyzed the data. L.B. lead the writing with contributions from S.S., T.J.N., Z.N. and J.R.-C.; Á.K. did the proofreading. All the authors have approved the manuscript.

**Funding:** The research was supported by the National Research, Development and Innovation Office (NKFIH; 108755). L.B. was supported by the ÚNKP-17-3 New National Excellence Program of the Ministry of Human Capacities. The project was supported by the European Union, co-financed by the European Regional Development Fund (EFOP-3.6.1-16-2016-00022).

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

#### **References**


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