Comparing the Runoff Decompositions of Small Experimental Catchments: End-Member Mixing Analysis (EMMA) vs. Hydrological Modelling

Abstract

This study is focused on the comparison of streamflow composition simulated with three well-known rainfall–runoff (RR) models (ECOMAG, HBV, SWAT) against hydrograph decomposition evaluated with End-Member Mixing Analysis (EMMA). In situ observations at two small mountain testbed catchments located in the south of Pacific Russia are used. All applied RR models and EMMA analysis demonstrate that two neighboring catchments disagree significantly on the mutual dynamics of the runoff sources. The RR models' benchmark test is based on proximity to EMMA hydrograph composition. Different aggregation intervals (season, month, and pentad) were applied to find a reasonable generalization period ensuring the clarity of results. ECOMAG is most conformable to EMMA outcome; HBV reflects flood events well enough; SWAT exhibits distinctive behavior compared to the other models. It is shown that, along with standard efficiency criteria of simulated and observed runoff proximity, EMMA analysis might provide useful auxiliary information for the validation of modelling results.

1. Introduction

Understanding runoff composition is a fundamental problem in hydrological science [1]. Most of the processes in hydrological systems take place underground and only a limited range of measurement techniques are available; the observation data need to be extrapolated in both space and time to predict the response of any catchment to a given rainfall event with an acceptable uncertainty [2]. Contribution of the different water masses to total runoff is time-dependent owing to the temporal variability of the flow response processes. In spite of the many publications on this subject, it is still a challenging task to determine how a catchment generates runoff as a whole system.

There is an undeniable consensus in the hydrological community that research on experimental catchments plays a key role in the understanding of the hydrological processes and provides essential outdoor laboratories for validation flow pathways and runoff generation mechanisms [3,4,5]. Small catchments are less inertial (more reactive to rainfall events) and offer potential noise reduction in the data compared with macro-scale data. Moreover, small catchments allow one to test hypotheses about runoff generation mechanisms due to relatively homogenous landscape properties. Alternatively, the greater catchment area complicates the interpretation of observational data and subsequent modelling results.

Sophisticated hydrological models are powerful tools to test hypotheses about catchment function and deal with uncertainties throughout the observation–conceptualization–modelling sequence [6]. Plausibly simulated streamflow can be obtained from models with distinct structures, providing the different runoff compositions, that is widely known as the equifinality effect. To narrow down the equifinality, it is necessary to evaluate the simulation results independently, that is, how realistically the model reflects the natural behavior of a hydrological system. However, the selection of hydrological models and the model structure uncertainties are still not fully understood [7,8,9,10,11]. Aggregation of outcome from different models or evolving model structure along axis of complexity using a “top-down“ strategy may lead to parsimonious model that provides useful insights into catchment behavior [10,12,13,14,15].

Catchment hydrology is still very much an empirical science [16,17], and it is common that some parts of a conceptual model may be more rigorously based on physical theory than others. Empirical approaches are still the basis or part of the numerous well-known rainfall–runoff (RR) models. The main reason to rely on this empirical knowledge is the scale-dependency of hydrological response unit (HRU)-based models, where they are more appropriate than any small-scale physical laws [18,19,20,21]. The model parameters are seldom directly measurable and are inferred from calibration. They are generally model-specific and represent average behavior in terms of both spatial and temporal variability [22]. Model structure uncertainty is as impactive on modelling outcome as parameter uncertainty, if mathematical formulation is not based on fundamental laws but adopts a set of lumped functions relating impact to the response [18].

All of the above creates a desire for alternative means of model evaluation in terms of how well a model captures the partitioning, storage, and release of water by a catchment [23]. End-member mixing analysis (EMMA) can provide estimates of the relative contributions of direct (surface) flow, deep groundwater flow, and soil percolation flow to total catchment runoff in the context of where these end-members can be easily separated and sampled [24,25,26,27]. Numerous researchers have indicated the importance of water chemistry for the transit time estimation and recognized the utility of tracers as an additional independent measure for model evaluation [27,28,29,30,31,32]. Despite gradual improvements in the reliability and economics of the field and laboratory methods, which provide necessary fine-time resolution data series at the present time [33,34], widespread application of EMMA is limited by the lack of suitable geochemical datasets.

This study is focused on a comparison of streamflow composition simulated with three well-known RR models (ECOMAG, HBV, SWAT) against the hydrograph decomposition obtained from EMMA. The main objective is to choose the RR model that best complies with EMMA in terms of the hydrograph separation that can be linked to a more accurate representation of real runoff generation processes. In general, it relates to the problem of harmonizing various approaches based on the solution of direct and inverse tasks of modelling. In the ideal case, the results of solving direct and inverse tasks should be the same or match closely.

The paper is laid out as follows: the next section describes the case study catchments, field observation details and all used models (EMMA, ECOMAG, SWAT, and HBV). Results of EMMA and hydrological simulations are presented in Section 3. For clarity, the results of hydrological model calibration are placed in the Appendix A, Appendix B and Appendix C. All results are generalized and discussed in Section 4; and finally, a summary and several concluding remarks are provided.

2. Description of Study Objects, Data, and Measurements Methods

The studied territory is the Pacific Russia boreal forests and is influenced by the East Asian Monsoon. Here, runoff modelling and hydrograph separation were performed for two small catchments (Elovy (3.5 km²) and Medvezhy (7.6 km²) creeks) that belong to the Upper-Ussuri Biocenological Experimental Station (45 km², 44°02′ N, 134°11′ E). The considered area is characterized by mid-mountainous relief with moderately steep (

Table 1. Catchments main characteristics.

Characteristics	Elovy	Medvezhy
Area, km2	3.5	7.6
Avg. height, m	722	704
Max. height, m	962	869
Avg. slope, %	13.5	13.8
Max. slope, %	28.7	31.5
Mean annual precipitation, mm 1	780	830
Mean annual temperature, °C 1	3.0	3.2
Avg. discharge, mm day−1 1	0.65	0.75

Note(s): ¹ based on observations for the period: 01.01.11–31.12.14 for Elovy creek and the period 01.01.14–31.12.17 for Medvezhy creek.

), locally very steep, hillslopesof up to 30% [3,35]. The average altitudes are 500–700 m a.s.l, and maximal values reach 1100 m a.s.l. (Figure 1).

Air temperature dynamics is characterized by high variability from year to year; average, and absolute minimal and maximal values are +0.7 °C, −45 °C (in January) +38 °C (July–August), respectively. The annual average precipitation amount is 700–800 mm, up to 80% of which occurs during the warm period (from late April to October) as rain. Unstable intra-annual and long-term precipitation regimes of the territory define runoff formation conditions of the investigated catchments. Resources of shallow water are not significant in terms of local water balance due to fractured rocks. Small-river runoff can be vanishingly low during winter or summer drought periods and reach 30–50 mm d⁻¹ during high-intensity rainfall caused by tropical cyclone–typhoon activity. The range of maximal daily heavy rain is 100–200 mm, and frequency of such events is assessed as one time per 5–6 years. The stable snowpack usually occurs in November, snow-cover depth reaches 0.5–1 m by the end of March, and the common range of the snow water equivalent is 100–200 mm.

Field observations were made from March to October in 2012–2018. Two weather stations (WS-GP1, Delta-T, Cambridge, UK) were mounted at 650 (Elovy catchment) and 750 m a.s.l. (Medvezhy catchment) to record air temperature and humidity, wind characteristics, and precipitation volume (see Table 1). Stream outlets were fitted with hydrostatic water-level loggers (Solinst 3001, Georgetown, ON, Canada). Observations temporal resolution was 15 min. Measurements of water discharges were run once in 2–4 days during low water periods, daily during average moisture conditions, and twice a day during flood periods. For discharge assessments, flow velocity was measured using magnetic inductive flow sensor (SEBA FlowSens, Kaufbeuren, Germany); cross-section geometry of riverbeds in permanent sites was measured simultaneously. To obtain daily discharge for the whole observation period, a water-level-discharge rating curve Q = f(H) was applied.

Hydrological properties of soil were obtained from 14 soil profiles (4 are located in the Medvezhy and 10 in the Elovy) [36]. The soil cover of the considered territory is dominated by the Cambisol group [37]. Brown forest soils here are characterized by a high rock fragment content (greater than 80% of profile volume) which leads to high infiltration rates (>30 mm/h).

Rainwater samples were collected to analyze chemical composition by using pyramidal funnels that were connected to polypropylene collectors. Soil water samples were taken with the vacuum lysimeters (DIK-8392, Daiki Rika Kogyo Co. Ltd., Kusatsu, Japan) installed at a depth of 0.35–0.80 m. Streamwater samples were collected simultaneously with discharge measurements. Streamflow temperature, pH, specific conductance, and total dissolved solids (TDS) were measured in situ using a handheld multiparameter meter (Yellow Springs Instrument, YSI Professional Plus, Yellow Springs, OH, USA).

The HCO₃ concentration was taken to be equal to total alkalinity, and estimated in the field laboratory after 2–4 h of sampling by the potentiometric titration method in unfiltered water samples using the standard technique. Water samples for chemical analysis were filtered at the field laboratory through the Durapore filter (Millipore, Burlington, MA, USA) with a pore size of 0.45 μm. To determine the content of cations, about 10 mL of each filtered aliquot was acidified with purified nitric acid and stored in polypropylene bottles at room temperature prior to analysis. Then, filtration samples of water for DOC analysis were stored in 30 mL glass bottles at a temperature of 4–8 °C [38,39,40].

Laboratory analysis of samples was carried out at the Geochemistry Laboratory of the Pacific Institute of Geography of the Far Eastern Branch of the Russian Academy of Sciences. The content of cations Ca, Mg, K, Na was determined by Atomic Absorption Spectroscopic method using a Shimadzu AA 6800 spectrometer. Cl, NO₃, SO₄ were determined using a Shimadzu LC 10Avp liquid chromatograph. Dissolved Si was determined by spectrophotometry method using blue complex with ammonium molybdate, and DOC by a TOC analyzer—Shimadzu TOC-VCPN.

Data quality control was performed to discard obvious errors associated with the operation of loggers. According to that, the next simulation periods were suggested: from 2012 to 2014 for Elovy creek, and from 2015 to 2017 for Medvezhy creek. The nearest WMO weather station 31939 (Chuguevka) is located 35 km northwest of the observation sites, so its data were used to fill up observational gaps and the cold part of the years. WRF-ARW model output was applied to derive daily solar radiation interpolated from closest computational grid nodes.

In this study, three RR models (SWAT, HBV and ECOMAG) represent a wide range of conceptual diversity of approaches used to solve the direct task of runoff modelling. HBV is lumped storage-based, and the other two are considered as spatially distributed (HRU-based) models. They differ in the structure and methods of runoff generation mechanism parameterization, and physical and mathematical foundation (see

Table 2. Hydrological models main characteristics.

Characteristics	ECOMAG	SWAT	HBV
Spatial discretization	HRU	HRU	Lumped
Temporal discretization	Daily	Daily	Daily
Number of calibrated parameters	11	12	10
Snow component basis	Degree-day	Degree-day	Degree-day
Potential evapotranspiration	Dalton method	Penman–Monteith	Penman–Monteith
Actual evapotranspiration	Linear reduction in PET by soil storage content	Reduction in PET by soil water content	Linear reduction in PET by soil storage content
Surface flow	Kinematic wave	SCS curve number	Linear storage
Soil flow	Darcy’s law	Kinematic storage model	Linear storage
Groundwater flow	Darcy’s law	Linear storage	Linear storage
Routing method	Kinematic wave	Variable travel time	Triangular weighted

).

The end-member mixing analysis (EMMA) framework is an example of solving the inverse task of modelling using methods of tracer hydrology. In contrast to imitation RR models, EMMA can be called an identification model, providing the interpretation and quantification of runoff generation sources. It is substantial: the discharges of different runoff sources are determined from direct observations of the water and tracer volumes using the mass balance equations system and some empirical relations. Overall, this could be interpreted as a sophisticated way to identify runoff sources.

All listed models are well known in the hydrological community and have previously been applied in the region [36,39,40,41,42,43]. For a description of the models’ structure, variables, and parameters, the original terms from the base publications are used. Detailed information can be obtained from references. A description of the models’ parameterization is provided in the Appendix A, Appendix B and Appendix C.

2.1. EMMA Model

End-Member Mixing Analysis (EMMA) is a commonly applied method to identify and quantify dominant runoff sources. This method is based on the assumption that waters, flowing through specific compartments of catchment, acquire its unique chemical signatures. Thus, some solutes can be used as natural tracers to determine the contributions of each water source (labelled as the “end-member”, EM) to total runoff. This hydrograph separation technique is based on the mass balances for the tracers and water, and assumes that: tracers are conservative; chemical signature of water contributions from various sources is constant and unique to make it distinguishable from others; mass conservation law applies to both the water quantity and solute quantity including mixing of water masses without their essential transformation [44,45].

Water runoff source identification and estimation of its contributions to streamflow (proportions) were performed using EMMA in combination with the diagnostic tools of mixing models. The algorithms for multidimensional statistical analysis used are implemented by the software library of extensions for Excel-Chemometrics [46]. According to the procedure described in Hooper [47], conservative tracers and the number of EMs that contribute to streamflow were identified from streamflow chemical data without using information about EMs. First, bivariate scatter plots of solutes against each other were analyzed for the linearity of mixing. Solutes that demonstrated collinear structure in the bivariate plots were considered as acceptable tracers. These solutes/tracers are variables of measured chemical data matrix G, n × p, where n is the number of samples, and p is the number of tracers. They are standardized using the mean and standard deviation of each solute. Standardized matrix G* was projected into a new Euclidean m-dimensional space U (m < p), obtaining matrix Ĝ* (n × m). This projection is made by principal component analysis (PCA) using eigenvectors S (p × p) extracted from the correlation matrix G*^T G* according to the equation

Ĝ* = G* × S^T × (S × S^T)⁻¹ × S

(1)

Conservative tracers and the number of EMs were determined by examining the relation of residuals E (the difference between Ĝ and G values, where Ĝ is destandardized Ĝ*) against the measured values for each solute. The dimension m of space U was increased until there was no structure (no correlation on the bivariate plots) to the residuals. Those solutes that exhibit a random pattern in the residuals, controlled by negligible correlation in the relations mentioned, are considered as conservative. Their number is denoted as k (k < p). The number of EMs is determined by m + 1.

Then EMMA was used with conservative tracers to identify EMs and to quantify the proportions of EMs in the streamflow following the procedure described in [45]. Streamflow samples chemical data were orthogonally projected using the truncated matrix of eigenvectors S₁ (k × m) extracted from the matrix G₁, which includes only k conservative tracers. Tracer concentrations in EMs were standardized using the mean and standard deviation of each tracer in streamflow samples and projected using the same eigenvectors S₁. Mixing diagrams were plotted in the U-space to screen EMs which contribute to the streamflow. In our study, the first two or three principal components (PCs), or the PCA scores, were used to determine the proportions of the EMs in the runoff composition. Thus, the PCs are considered as the tracers for the mixing model—two tracers for a three-source mixing model or three tracers for a four-sources mixing model.

Determining the proportions of sources in the streamflow based on the mixing model is performed in the usual way, by solving a system of three or four mass balance equations for water and tracers, as described in [45]. Model validation was performed by multiplying the hydrograph separation results (proportions) on the measured solute concentrations in the EMs to reproduce the chemical data of streamflow for conservative solutes. The high correlation between the series of modelled and measured tracer concentrations in streamwater indicates the adequacy of the obtained mixing model.

Hooper [47] suggested that assumptions of linearity of mixing and conservative behavior of tracers can be evaluated using bivariate scatter plots and residuals derived from the selected model. Bivariate scatter plots should be developed for all potential combinations of available solutes. A collinear structure in the bivariate plots could be used to infer conservative behavior. These data were used to construct a correlation matrix followed by principal components analysis (PCA) to extract eigenvectors and eigenvalues. The eigenvectors form the basis for a new Euclidean space, U-space.

2.2. ECOMAG Model

ECOMAG (ECOlogical Model for Applied Geophysics) is the process-based, semi-distributed RR model [48]. The spatial structure of the ECOMAG model splits watershed into sub-basins based on topography, river network structure, soil and vegetation type, land use, and variability of climate characteristics. The main ECOMAG model equations were adopted from the full spatially distributed model [49] by neglecting secondary terms and spatial aggregation at subbasin scale. Daily resolution time series of precipitation, air temperature, and air relative humidity are used as inputs. Computation of river basin hydrological response is described by two main phases: calculating the effective precipitation for each sub-basin and then routing it through the river network. Runoff from sub-basin is calculated as the sum of three components: Horton overland (surface) flow (when rainfall has exceeded infiltration capacity and depression storage capacity), soil flow (sum of runoff from horizon A and B) and groundwater outflow (baseflow and infiltration from horizon B). In warm periods, precipitation is partially infiltrated and moves along the hillslopes (over impermeable surfaces) as interflow. Excess water produces surface flow and moves downslope towards the drainage network. The rest of the water that has not been drained into rivers as lateral or surface flow can be evaporated or percolated into deep aquifers. Within cold and mid-season periods, the model describes snowpack evolution and soil freezing–thawing cycle.

2.3. SWAT Model

The Soil and Water Assessment Tool (SWAT) is a semi-distributed hydro-ecological model [50,51]. This model uses the net of hydrological response units (HRUs) as the base elements to model main hydrological cycle processes—infiltration, evaporation, runoff generation, and soil hydrothermic and snow-cover dynamics. The daily volume of precipitation excess and surface runoff is calculated by the empirical SCS CN method. Kinematic wave approximation is used for channel routing. Daily rainfall, maximum and minimum air temperature, solar radiation, relative air humidity and wind speed are the inputs for the model. The SWAT uses air temperature to define input precipitation as rain or snow. Some of the precipitation, intercepted by canopy, can be evaporated. Then, the simulation is performed as two successive land and routing phases. The snow thawing is linearly related to snow cover depth and fraction of watershed area and air temperature. At the top of the soil profile, precipitated water can either flow overland or infiltrate into underlying strata if their temperature is positive, the field capacity of the upper layer is exceeded, and the underlying soil layer is not saturated. Water may be partially evaporated and turned into lateral or groundwater runoff components inside the soil profile. The dynamic of the lateral flow runoff component is calculated using the saturated conductivity of each soil layer. Percolated through a soil column, water is proportionally divided between unconfined and deep aquifers. Potential evaporation is used to calculate the water exchange between groundwater and the bottom of the soil profile.

2.4. HBV Model

The conceptual HBV model was developed by the Swedish Meteorological and Hydrological Institute (SMHI) [52]. The model consists of three main modules: snowmelt and snow accumulation, soil moisture and effective precipitation routine, and runoff transformation to catchment ’s outlet. The measured precipitation is supposed to be snow if the air temperature is lower than threshold temperature TT; otherwise, precipitation appears as rain. The simulated snowpack volume can be adjusted with correction factor SCF. The degree-day method is used for snowmelt calculation. Groundwater recharge and actual evaporation are simulated as functions of actual soil storage. The amount of water available for runoff generation is calculated as the ratio of actual soil moisture SM to field capacity FC and the power coefficient Beta-(SM/FC)^Beta. The model first recharges the upper storage and then the lower storage using percolation parameter PERC. Surface flow appears when the upper storage capacity exceeds a certain threshold. Potential evapotranspiration (PET) is input to the model, which can be tuned by the Cet parameter. Actual evaporation is equal to PET when SM/FC is higher than the LP parameter; otherwise, a linear reduction is used. The runoff is calculated as the sum of three linear outflows: surface flow Q0, interflow Q1, and baseflow Q2 with three correspondent recession coefficients K0, K1, and K2. A transformation function with the triangular weighting parameter MAXBAS [53] is used for smoothing the total runoff to obtain the discharge at the outlet.

3. Results

3.1. EMMA Results

The EMMA procedure was performed at the Medvezhy catchment for the period 2015–2017, and at the Elovy catchment for the period 2012–2014. Bivariate solute plots were constructed for all pairs of chemical indicators using streamwater samples for the whole period of observations. Diagrams depicting collinearity and linear trends were used to define tracers with conservative behavior. Stronger linear trends between pair solutes were observed at the Medvezhy catchment, and weaker linear trends were observed at the Elovy catchment, R² < 0.5 (see

Table 3. Pairwise correlations of the concentrations of solutes, ranked by R² (n is the number of samples).

R2	Medvezhy Catchment, n = 69	Elovy Catchment, n = 126
>0.71	TDS—HCO3, SO4—Mg, HCO3—Mg, TDS—Mg	–
0.7–0.61	DOC—SO4, HCO3—SO4, SO4—Ca, HCO3—Ca, TDS—Ca, Ca—Mg	–
0.6–0.51	TDS—SO4, DOC—Mg	–
0.5–0.41	–	HCO3—Mg, HCO3—TDS, SO4—NO3
0.4–0.3	–	NO3—HCO3, NO3—Na, HCO3—Na, Na—TDS

).

For the Medvezhy catchment, PCA was performed for the matrix of chemical indicators including six solutes: TDS, DOC, HCO₃, SO₄, Ca, and Mg. The first two PCs explain more than 92% of the total variance of these data (see

Table 4. The percent of variance explained by PC (U) for individual catchments.

Catchment	Tracers	U1 *	U2	U3	U4	U5	U6
Medvezhy	TDS, DOC, HCO3, SO4, Ca, Mg	83.6	8.6 (92.2)	4.4 (96.7)	2.1 (98.7)	0.9 (99.6)	0.4 (100)
Elovy	HCO3, Mg, NO3, Na	64.8	20.2 (85.0)	9.7 (94.7)	5.3 (100)	-	-

Note(s): * Values in parentheses—the cumulative percent of total variance explained by set of PCs from the first to the current.

). Analysis of PCA-model residuals against measured values for individual solutes at two-dimensional mixing subspaces is shown in Figure 2a. The random pattern for each solute indicates that these six tracers could be considered as conservative and can be used for adopting the mixing model with three EMs.

For the Elovy catchment, a four-source mixing model is required to successfully represent the runoff generation in the Elovy Creek [39,54,55]. The fact that HCO3, Mg, NO3, Na tracers captured the 95% of the total variance is explained by the first three principal components. The hypothesis about conservatism of tracers (stable concentrations of tracers along their pathways) is confirmed by the randomness structures of residuals shown in Figure 2b.

Here, the strict randomness of the residuals is not always confirmed. Several known tests for randomness do not seem quite suitable because of inconvenient properties of data series—non-Gaussianity, uneven sampling, and discontinuous intra-series connectivity. This uncontrollably reduces the reliability of the test assessments. In this case, it is advisable to use only the simplest criteria and use their conclusions only as one of the arguments. We used as a test the critical values of the R² assessments by sample (with its true value zero), calculated based on the well-known Fisher’s normalizing z-transform at a significance of 99% [56]. For samples of 69 and 126 members, the critical values of R² are approximately 0.101 for the Medvezhy and 0.054 for the Elovy, respectively (compared with Figure 2).

In the case when the residuals are not completely random, but with an insignificant value of R², the inclusion of the tracer can be justified by independent considerations and analysis of the simulation results as a whole. For example, R² for Ca (see Figure 2a) is essentially non-random. However, there is a group of three outliers in this residual plot, related to samples taken during one week in May 2015. The exclusion of this group from the sample reduces the value of R² by almost half, which serves as an argument for including this tracer in the model. The use of not completely conservative tracers decreases the accuracy and reliability of estimates in the mixing model. This fact was taken into account when discussing the conclusions below. The example above probably indicates that the existing dataset does not sufficiently account for the seasonal dynamics; hence the mixing models are considered reliable only for the summer-autumn rain-flood season (from June to September for the studied territory).

Mixing diagrams were constructed for each creek by orthogonally projecting the conservative tracer’s matrix of river samples into U-space. Tracers that averaged concentrations of EMs (

Table 5. Chemical composition of end-members. n—number of samples, range of concentrations is in the numerator, average values are in the denominator.

Tracer, mg/L	Medvezhy Catchment			Elovy Catchment
	End-Members
	Rainwater (n = 36)	Groundwater (n = 15)	Soilwater (n = 12)	Rainwater (n = 36)	Groundwater (n = 5)	Soil Water 1 (n = 13)	Soil Water 2 (n = 7)
HCO3	<0.1–2.30	75.6–90.3	14.6–24.3	<0.1–2.30	14.2–19.2	10.4–18.0	6.71–9.76
	0.37	84.3	19.4	0.37	16.4	13.7	8.11
Mg	<0.02–0.13	1.66–5.18	3.25–5.36	<0.02–0.1	0.85–1.04	0.34–0.61	0.39–0.55
	0.05	3.47	4.56	0.05	0.98	0.4	0.48
TDS	0.05–14.9	121–163	30.5–100	-	-	-	-
	5.32	146	63.8
DOC	0.18–5.60	2.50–4.76	4.10–42.9	-	-	-	-
	1.77	3.56	18.1
SO4	<0.32–5.75	16.2–29.0	4.07–9.09	-	-	-	-
	1.42	23.9	6.71
Ca	<0.1–6.78	17.8–38.9	6.96–22.7	-	-	-	-
	0.65	27.8	15.4
NO3	-	-	-	<0.2–5.1	<0.2–0.72	<0.2–0.15	2.13–7.45
				1.17	0.39	0.05	4.66
Na	-	-	-	0.04–1.0	2.82–2.99	4.72–7.80	2.00–2.66
				0.14	2.91	5.79	2.41

) were also projected into the same U-space using eigenvectors. The two-dimensional mixing diagram for Medvezhy creek is presented in Figure 3a. River samples are located in the mixing subspace of tracers <U1, U2> and are bound by three sources. The mixing diagram for Elovy catchment represents the cloud of streamflow samples, which are located in the three-dimensional mixing subspace <U1, U2, U3> and determined by four sources (see Figure 3b). In both cases, the EM-based figures (triangle or tetrahedron) enclose most of the streamflow samples in the U-space. This geometrically confirms the possibility of representing each streamflow sample composition as a mixture of runoff sources, and generally verifies the mixing model.

The identification and interpretation of EMs for Medvezhy creek is unambiguous. Two EMs are related to ground and soil water (GW and SW1), and the third one to rainwater (RW). The GW is associated with streamflow samples taken during the low flow period and represents the baseflow component of the runoff. Samples of SW1 were taken from the upper organic horizon of soil and can be interpreted as lateral (or inter-soil) flow. RW reflects quick overland runoff that reached the drainage network almost without chemical transformation. The mixing subspace for Elovy catchment evolved to four EMs. The first three of them (RW, GW, and SW1) are the same as for Medvezhy, and the fourth, SW2, is associated with soil mineral source and sampled with a lysimeter at the bottom of the steep hillslope from fir-spruce forest soil with slow organic matter destruction.

The mixing models were finally verified by direct calculation of modelled tracer concentrations for every streamwater sample using EM tracer concentrations and estimated proportions of water runoff sources. The modelled and measured series of streamwater sample tracer concentration (Figure 4) are characterized by relatively high R² values (0.75–0.98).

The EMMA runoff components and measured discharges were used to build the relationships: GW = f (Q) and SW1 = f (Q) at Medvezhy creek; GW = f (Q), SW1 = f (Q) and SW1 + SW2 = f (Q) at Elovy (see Figure 5). The proportions of RW (for both streams) and SW2 (for Elovy), which are weakly related to measured discharges, are determined as the remainder from the water balance equation. These relationships were used to calculate the daily streamflow composition series using the measured hydrographs.

3.2. Results of Hydrological Simulations

All used models were calibrated independently. Hydrological simulations were made continuously (winter season included) with daily time step and a one-year warm-up period. Calibration was performed manually using a trial-and-error approach. The details of the models’ parameters values appear in the Appendix A, Appendix B and Appendix C. Observed and simulated hydrographs for typhoon-induced flood events of 2012 and 2016 are shown in Figure 6. Model efficiency was assessed using common goodness-of-fit measures (determination factor R², Nash and Sutcliffe efficiency (NSE) [57], and relative bias (BIAS, %) against the measured discharge at catchment outlets (see

Table 6. Goodness-of-fit characteristics for runoff simulations with different models.

Model	Medvezhy				Elovy
	Year	R2	NSE	BIAS, %	Year	R2	NSE	BIAS, %
ECOMAG	2015	0.92	0.91	−5	2012	0.93	0.91	−0.3
	2016	0.92	0.90	8	2013	0.87	0.80	13
	2017	0.91	0.87	−6	2014	0.84	0.83	4
	2015–2017	0.92	0.90	4	2012–2014	0.90	0.89	10
SWAT	2015	0.94	0.88	−20	2012	0.90	0.90	−4
	2016	0.89	0.85	6	2013	0.82	0.81	−1
	2017	0.69	0.67	−12	2014	0.93	0.88	−13
	2015–2017	0.90	0.86	−1	2012–2014	0.86	0.86	−6
HBV	2015	0.35	0.35	−8	2012	0.96	0.96	1
	2016	0.92	0.91	18	2013	0.83	0.83	−6
	2017	0.56	0.53	−3	2014	0.64	0.62	4
	2015–2017	0.91	0.91	12	2012–2014	0.88	0.88	−0.3

). NSE is categorized as “very good” when its value > 0.75 and “unsatisfactory” when its value < 0.5, interim ranges (0.75 > NSE > 0.65 and 0.65 > NSE > 0.5) are defined as “good” and “satisfactory”, respectively. According to BIAS values, simulation results are assumed to be unacceptable if |BIAS| > 25%, “satisfactory” if |BIAS| > 15% and < 25%, “good” for if |BIAS| < 15% < and > 10%, and “very good” if |BIAS| < 10% [58]. According to these criteria, more complex and sophisticated models provide better results for both catchments, and ECOMAG outperforms SWAT in most cases. HBV provides good results for the period of high flow but demonstrates poor performance for low flow (0.01–1.0 mm d⁻¹) periods.

4. Discussion of Hydrograph Separations

The dynamics of runoff composition for the years with significant flood events is presented in Figure 7. All RR models demonstrated that two neighboring small catchments significantly differ in the mutual dynamics of the runoff sources. Geochemical analysis confirmed the differences in runoff generation processes at the both studied catchments. For the Medvezhy catchment, the mixing model provided a typical three-source stable structure and reliable estimates of the source proportions. The proportion of soil runoff for this catchment is relatively small. Conversely, the soil runoff proportion for the Elovy catchment is more significant against surface and groundwater. The younger geology and higher elevation range within this smaller and narrower catchment, with at least two apparent soil–vegetation belts, led to applying a four-source (with two different soil water sources) EMMA model.

From the resulting simulated series, data on runoff constituents were extracted, which directly depend on the model applied. These data were compared with runoff sources obtained from EMMA. Here, we would like to highlight the two main aspects concerning the correctness of further comparison of runoff composition provided by different models.

First, it is difficult to reconcile the interpretations of modelled runoff constituents. The hydrological models derive runoff sources from hydro-physical and hydraulic characteristics of each proposed flow path. Their results are conditioned by the different basic concepts used, always more or less speculative. On the other hand, the tracer separation model is based on the interpretation of measured hydrochemical data related to water residence time, geochemical features of rocks, and interaction intensity with soil matrix. The runoff sources in different models, denoted by the same terms, may not have the same sense. The problem of compatibility and reconciling of different interpretations of the river flow structure is one of the most complex problems of runoff theory and is beyond the scope of this study.

Second, the comparison may be complicated by relatively low accuracy and instability in daily runoff proportions, “measured” by EMMA. This is due to both the well-known limitations of EMMA (incomplete conservatism of tracers and spatial variability of sources) and the limited number of samples leading to significant dispersion in statistical relationships used to calculate daily series of the runoff constituents. Concerning this latter problem, data aggregation in longer time intervals should allow one to mutually compensate for the random errors and to smooth the instability in daily runoff fractions. Therefore, the range of data aggregation intervals (season, month, and pentad) are applied to find a reasonable data generalization period, ensuring the clarity of the results.

The further benchmark analysis is based on the ranking of proximity of runoff constituents simulated by RR models to hydrograph decomposition by EMMA, regarded as an identification model. Each RR model used (ECOMAG, HBV, SWAT) provides three runoff fractions related to surface, soil, and groundwater. The EMMA mixing model separated the Medvezhy catchment runoff to three sources that correspond to the runoff structure of hydrological models. For the Elovy catchment, two EMMA soil flow sources are summarized to match the structure of the three runoff components’ structure in other models.

As the first step, the most generalized seasonal (June–September) data aggregation (see

Table 7. Seasonal runoff constituents obtained from rainfall–runoff models and hydrograph separation by EMMA.

Source	Medvezhy				Elovy 1
	SWAT	HBV	ECOMAG	EMMA	SWAT	HBV	ECOMAG	EMMA
	2015				2012
Surface, %	0.0	0.0	43.5	25.6	0.0	22.1	7.7	15.1
Soil, %	97.8	2.8	29.9	14.6	35.3	72.3	89.3	75.5
Ground, %	2.2	97.2	26.6	59.8	64.7	5.6	3.0	9.4
Total, mm	26.9	43.9	34.6	37.2	64.6	74.1	78	74.1
	2016				2013
Surface, %	57	7.8	32.1	34.4	0.0	5.0	3.5	24.4
Soil, %	33.3	43.7	63.0	33	36.2	85.9	94.2	66.3
Ground, %	9.7	48.5	4.9	32.6	63.8	9.1	2.3	9.3
Total, mm	211.3	233.4	213.3	195.2	62.2	76.6	100.3	85.6
	2017				2014
Surface, %	0.7	0.0	58.3	27.3	0.1	1.0	3.8	23.8
Soil, %	99.3	13.3	22.9	18.3	45.4	86.2	92.8	63.2
Ground, %	0	86.7	18.8	54.4	54.5	12.8	3.3	12.9
Total, mm	35	38.1	37	38.7	49.4	73.7	71.0	66.6
	2015–2017				2012–2013
Surface, %	44.1	5.7	35.8	32.1	0.0	9.3	4.9	21.2
Soil, %	47.8	34.4	54.8	28.4	38.4	81.5	92.3	68.4
Ground, %	8.1	59.9	9.4	39.5	61.6	9.2	2.8	10.4
Total, mm	273.1	315.5	284.9	271.2	176.2	224.3	249.3	226.3

Note(s): ¹ Soil flow proportion for Elovy catchment obtained by EMMA is the sum of the SW1 (soil organic) and SW2 (soil mineral) parts; see Section 3 for details.

) was examined. The total seasonal runoff volume for the Medvezhy catchment is reproduced more accurately: by ECOMAG in 2015 and by HBV in 2017; in 2016, SWAT and ECOMAG demonstrate almost identical proximity. The seasonal runoff composition, provided by EMMA, was captured best using ECOMAG for the dry seasons of 2015 and 2017. However, during the high-water season of 2016, it was impossible to choose the best hydrograph separation (as the closest to reference EMMA) among the three simulated variants discussed.

In most cases of the Elovy catchment runoff separation (see Table 7), HBV seems slightly better than ECOMAG, and it was the most successful in the estimation of the groundwater constituent. SWAT for both catchments provided an essentially different runoff composition from EMMA and greatly underestimated the Elovy catchment total seasonal runoff volume. Thus, seasonal data aggregation is not optimal for an unambiguous selection of the best simulation model, at least for the available data.

Next, monthly data aggregations were used to rank the RR models according to a better approximation to the reference EMMA. The monthly aggregation of the three runoff simulation results was compared with 12 observed monthly runoff volumes (four months in each of the three years) and their 36 source proportions estimated by EMMA (see Figure 8). For each case, the simulation providing the closest approximation was fixed, and then the total number of the best results was calculated.

For the Medvezhy catchment, the closest proximity of the ECOMAG and EMMA runoff composition is more obvious (see Figure 8). Comparison of monthly runoff volumes showed the approximate parity between HBV and ECOMAG-7 versus five cases (SWAT loses all the cases). ECOMAG runoff composition is close to EMMA results in 22 out of 36 cases; it is followed by HBV with 11 cases, while SWAT is successful only in 3 cases.

Conversely, for the Elovy catchment, the runoff compositions of the HBV and ECOMAG models are close to the EMMA hydrograph separation (see Figure 8). For runoff volumes of the Elovy catchment, both ECOMAG and HBV score five cases, and SWATthe remaining two cases. The Elovy monthly runoff source proportions provided by HBV and ECOMAG are close to EMMA estimates in 16 cases, and SWAT—in 4 cases. An advantage of the HBV model clearly appears in the estimation of the groundwater fraction in streamflow. This is also an advantage of the ECOMAG model—in the estimation of the surface water fraction. The SWAT model hydrograph decomposition based on monthly analysis appear to be uncompetitive.

Based on results of a five-day data aggregation comparison (see Figure 9), the advantage of the ECOMAG model becomes more pronounced for the Medvezhy catchment. The ECOMAG runoff source dynamics is very similar to the EMMA decomposition. However, for the outstanding 2016 flood event, HBV results are closer to EMMA, while the ECOMAG provided significantly different runoff composition. For the Elovy catchment, it is hard to choose with confidence the best approximation to EMMA decomposition from the results of the ECOMAG and HBV models. It should be noted along the three years (2012–2014) forming the observation series on the Elovy catchment are similar in terms of volume and regime of runoff, in contrast to 2015–2017, included in the observation series on the Medvezhy catchment. An extreme event similar to the flood of 2016 was not observed at the Elovy catchment.

In accordance with comparison results, the ECOMAG model yielded the best agreement with EMMA, HBV was second, and the SWAT model (despite outperforming HBV by NSE efficiency somewhat) was last. The competitive advantage of the ECOMAG model is provided by the proportion comparability of all three runoff sources and the mutual dynamics of their proportions as well. However, while surface runoff estimated by ECOMAG is in good agreement with EMMA, groundwater source estimations differ markedly. In some cases, the HBV model is able to reproduce the EMMA ground runoff quite accurately, but rather the overall pattern over relatively long (seasonal) intervals, not its detailed behavior. At all times, SWAT showed more difference from EMMA than other models.

EMMA clearly explained the differences in hydrological behavior between the catchments due to the need to use four end-members for the Elovy catchment instead of the usual three for the Medvezhy catchment. The need to use a four-component mixing model suggests that the soil interflow of Elovy catchment is more complex in its organization compared to that ofthe Medvezhy catchment.

The simplified lumped models (e.g., HBV) provide the classical storage-based composition of the resulting hydrograph. More sophisticated models, by means of the hydrological characteristics of the soil horizons, allow one to separate the interflow by means of soil horizon hydrological constants (WP, AWC, FC), and vertical and horizontal filtration within one or more layers.

All the above-noted inferences about the representative features of runoff generation in various models were formulated only preliminarily based on the limited data used. The difference in performance can also be partly attributed to the subjective choice of the final model parameterization. It should be concluded that the monthly and synoptic (five/ten days) data aggregation are the scales at which the seasonal runoff dynamics is well suited to clearly compare the RR models with representation of flood events in aggregated view.

5. Conclusions

This article demonstrates the results of comparing the simulated runoff hydrographs and runoff constituents obtained by calibration of three different RR models (SWAT, HBV, ECOMAG), with observed runoff hydrographs and the runoff sources’ dynamics provided by EMMA. For the study, we used in situ observations in two typical small experimental catchments of the mountain-taiga region in the south of Pacific Russia. Both the RR and tracer models provide unequivocal evidence that two neighboring small catchments significantly differ in the structure of runoff generation processes related to differences in the geological and landscape structure. It is shown that compliance with EMMA results may provide useful auxiliary information for validation of hydrological modelling results besides the use of standard efficiency criteria that reflect the proximity of simulated and observed runoff values.

Summarizing the above, we can give a preliminary answer to the question that motivated this study. Comparison of runoff composition obtained with direct and inverse modelling can rank the RR models in terms of physical description completeness of the runoff generation processes. At the same time, the problem of reaching a detailed model-to-model correspondence of runoff composition dynamics remains an urgent challenge. This requires development of unified classifications of water mass types within a river basin and knowledge of runoff generation mechanisms taking into account regional specifics. All conclusions are preliminarily based on limited observation data from two analyzed objects. The main perspectives of the research require progress towards development in methodology, terminological, and conceptual compatibility of interpretations of streamflow composition behind various approaches.