The specification of the model consists of transforming the researcher’s perceptions into the formal configuration of a structural equation model. This transformation is the most important and complex task in SEM because each of the following steps is grounded on the premise that the model designed is properly specified and that only a correctly specified model can properly test the researcher’s hypotheses (Ntoumanis and Myers, 2016).

According to the literature review presented in the Introduction, fluvial disturbance directly impacts riparian vegetation in two major ways, mrphd and physiological stress. A mrphd means surface flow-derived processes that cause physical vegetation damage, sediment burial or uprooting, and entrainment (e.g., Friedman and Auble, 1999; Edmaier et al., 2011; Bendix and Stella, 2013). Physiological stress is a consequence of the grndh and impacts riparian vegetation as water stress caused by lowering of the level of groundwater table during a drought or by anoxia during flood periods (e.g., Stromberg et al., 1996; Shafroth et al., 2000; Cooper et al., 2003; Baird et al., 2005; Lowry et al., 2011). The current patch mosaic (i.e., the location and shape of the patches) of riparian vegetation is therefore the outcome of the existing conditions in the habitat, mainly as a consequence of the historical flow regime (Pettit et al., 2001). Consequently, we based ourselves on the literature review in the Introduction and our expert knowledge of riparian ecosystems to conceptualize two theoretical model constructs to address these research questions (Figure 1).

Morphodynamic disturbance and groundwater hydrology are the exogenous variables that depict the flow regime and have a direct effect on the endogenous variables patch location (phslc; Figure 1A) and patch shape (phssh; Figure 1B). These variables cannot be measured directly, so they are considered latent variables expressed by manifest variables that are measurable. Because both mrphd and grndh are controlled by the flow regime, they are expected to be correlated to some extent. The physiological effects of the duration of flooding on vegetation were not expressly accounted for by the models because they were not expected to play an important role in the riparian ecological succession in the Mediterranean watersheds considered in this study. The characteristic flashiness of the pluvial flow regime defines very short flood durations in small catchment areas, and furthermore, are restricted to the winter (Tockner et al., 2009), which is the dormant vegetation season, when floods have less of an effect (Crawford, 2003).

Finally, the model constructs were hypothesized to be parsimonious as possible, and the latent variable indicators were reduced by the maximum amount. In fact, single or only a small number of latent variable indicators are normal in the natural sciences due to the nature of the data (Grace, 2006). This strategy was adopted not only to avoid identification problems (Schumacker and Lomax, 2010) but also to make use of only one or two of the best indicators as recommended and often sufficient (Hayduk and Littvay, 2012), as well as for the sake of the simplicity that must be sought for any ecological model (Jackson et al., 2000).

The objective of model identification is to determine if a theoretical model construct enables the unique estimation of the requisite parameters from existing non-redundant information in the data. Consequently, to assure the identification of the proposed models, the factor loadings of the latent variables with only one indicator were set to 1 and the corresponding indicator error variance was set to 0 (Beaujean, 2014). Hence, 4 factor loadings exist for the specified models, 4 measurement error variances, 2 path coefficients, 1 correlation between the latent variables and 3 equation error variances – a total of 14 parameters that must be estimated. Each of the model constructs has only 8 free parameters (relationship coefficients to be estimated from the collected data) and 4 observed variables, implying that there are (8 × 4)/2 = 16 pieces of non-redundant information. Because only 14 parameters must be estimated, the models are considered overidentified and therefore the model identification is verified. Additionally, to have sufficient variability to estimate the model parameters, the sample size (N) should follow the rule of at least 20 cases for each free parameter that must be estimated (Jackson, 2003). Accordingly, in this situation N should be at least 160 observations. Notwithstanding, the usual minimum sample size for SEM studies is approximately 200 cases (Kline, 2011).

Four study sites were selected in natural conditions of riparian vegetation and flow regime (Figure 2). Upstream of the study sites, the main land uses in the watersheds are planted forests and natural shrublands, with very sparse villages and no noteworthy industry. In all cases, the flow regime was considered unregulated and typically Mediterranean, with a low winter flow interspersed by flash floods, and a very low and often intermittent summer flow (Bonada and Resh, 2013). The woody riparian species composition was similar in the study sites, comprising mainly willows (Salix spp.) and ashes (Fraxinus angustifolia). Notwithstanding, the four sites encompassed diverse fluvial geomorphologies, watershed features and river sizes (Table 1).

A complete survey of the topography and riparian vegetation was performed at each study site. The surveys were done on river reaches 300–500 m long (depending on the river width maintaining a ratio between length and main channel width of 10–20) and extended laterally to the area flooded by a 100-year recurrence interval (normally determining study site widths of approximately 70–110 m). This area was determined by modeling such a flood in each study site using the hydrodynamic model River2D. The topography was surveyed with an elevation detail of 20 cm using both total stations (Leica TPS400) and Global Navigation Satellite Systems (Leica 500 GPS, composed of two double-frequency-to-real-time SR 530 RTK antennas L1 and L2 AT 502). Riparian vegetation surveys were performed using a sub-meter handheld GPS (Ashtech, Mobile Mapper 100) to outline and georeference all the existing vegetation patches. Vegetation inventories in each patch sought to characterize succession phases by vegetation attributes and included phanerophyte species identification (Table 1). Then, each vegetation patch was classified by its succession phase, based on the vegetation type and the patch age. The succession phase classification methodology followed Rivaes et al. (2013) in which vegetation types were defined by indicator species and patch age by dendrochronological methods. A total of 220 vegetation patches in the aquatic, bank and floodplain zones were assessed. Four succession phases were found at all the study sites, namely, initial phase (IP), pioneer phase (PP), early successional woodland phase (ES), and established forest phase (EF). Altogether, those succession phases account for the existing ecological succession phases of riparian woodlands in those Mediterranean rivers (García-Arias et al., 2013). In detail, IP is characterized by areas with less than 50% vegetation cover and the absence of woody species. PP is typified by the recruitment areas of woody species and ES is characterized by the presence of well-established softwood pioneer individuals such as willows. EF is found in patches presenting hardwood species such as ashes, along with a well-defined understory stratum (Table 1).

The hydraulic parameters were obtained with the River2D model (Steffler et al., 2002) and were used to create the shear stress maps of the annual maximum instantaneous discharges and the mean annual water table elevations. The River2D is a hydrodynamic 2D model based on the depth averaged Saint Venant equations that computes the depth and the discharge intensities in the x-y directions. This tool was developed for application in natural rivers and features wet-dry area solution capabilities by combining surface flow and groundwater flow equations to compute the free elevation above and below the ground. River2D also incorporates a bed resistance model and a transverse shear model. In the former model, bed shear stresses are assumed to be related through the effective roughness height to the magnitude and direction of the depth-averaged velocity. The advantage of using the roughness height as the resistance parameter is that it remains constant over a wide range of depth. In the transverse shear model, the depth-averaged transverse turbulent shear stresses are modeled using a Boussinesq type eddy viscosity formulation. A complete description of the model is provided by Ghanem et al. (1996). Shear stress has been widely used as a fundamental proxy for soil erodibility, mrphd and drag imposed on vegetation by river flows. In natural channels, shear stress is considered to be balanced by three resistance components, namely, viscous drag on the ground surface on particles, pressure drag associated with large non-vegetal boundary roughness and drag on vegetal elements (Temple et al., 1987). The maps produced by River2D had a precision of a quarter of a square meter and were used to compute the MDi and the GWDi developed by Egger et al. (2014, 2016). Both are time and intensity weighted indexes (TIWI) that characterize the processes of grndh and mrphd of the historical flow regime in the past decade at each study site. These indexes are proxies for the long-term ecosystem processes of physiological stress and geomorphic-mechanical disturbance that effect riparian communities and provide a parameter that allows a dynamic analysis of riparian ecosystem patterns. Mean values of MDi and GWDi were calculated for the area covered by each of the vegetation patches recorded. The patch features were obtained using the ArcGis 9.2 (Environmental Systems Research Institute, 2010) and its Patch Analyst extension (Rempel et al., 2012). The relative positioning of the patches in the study sites were characterized by height (THAH) and distance (THAD) to thalweg. The patch shape metrics selected were the patch perimeter (PERIMTR) and the mean patch fractal dimension (MPFD). PERIMTR is an edge metric that provides a measurement of the dimensions and amount of edge created by each patch. This variable significantly affects many ecological phenomena and the analysis of spatial patterns in landscape ecological research (McGarigal and Marks, 1994). The MPFD is a shape metric that represents the geometric complexity of the patches and provides information regarding the formation and quality of the patches (Imre and Bogaert, 2004). Finally, the patch data were compiled and a data matrix was built and uploaded in the R environment (R Development Core Team, 2011) for subsequent data validation and treatment.

The variables were log transformed to remove bias and to bypass model under-identifiability imposed by disproportionate scales (McDonald and Ho, 2002). Outliers for which no justification was found were removed. Data were subsequently checked for compliance for the statistical assumptions for the SEM. Data normality was assessed with the Shapiro–Wilk normality test (Shapiro and Wilk, 1965), supported by an assessment of data sk and ku. Sk and ku were also used to assess the data normality to some extent, providing information about the distribution of the variables, i.e., the symmetry and the “peakedness” of the distribution. In this particular case, even after the data were transformed, some of the variables considered continued to have normality issues and did not pass the Shapiro–Wilk normality test. However, the magnitude of the sk and ku of the variables were always below 0.88 and 1.01, respectively. No consensus exists as yet about the maximum magnitude of sk and ku that does not undermine the reliability of the conclusions about the model fit and the parameter estimates (Finney and DiStefano, 2006) but the most conservative thresholds found in the literature are maximum absolute values of 1 and 1.5, respectively (e.g., Schumacker and Lomax, 2010; Kline, 2011). Furthermore, for large sample sizes (N ≥ 40) non-normality is not considered problematic (Pallant, 2002) and can be ignored (Altman and Bland, 1995). Consequently, the data were considered to be normal. The linearity between the variables was assessed using the Harvey-Collier and Rainbow tests for linearity, as well as Ramsey’s RESET test for functional form. The majority of the variable relationships passed all three tests and only two relationships were linear in at least one of the tests with a 99% confidence level. Consequently, the variables were considered to be linearly related. Independence was tested using Moran’s autocorrelation coefficient (Moran, 1950) and the Mantel test (Mantel, 1967) for spatial correlation. The data were considered independent, because both Moran’s coefficient and the Mantel test did not demonstrate any spatial correlation of the riparian succession phases between sites or at each site, with a confidence level of 99%. Homoscedasticity was assessed by a plot analysis of the residuals versus the predicted values. No substantial problems were found by a visual assessment of the plots and therefore, the data was considered homoscedastic. Multicollinearity was evaluated by calculating the variance inflation factor (VIF) for each variable. The VIF is a widely used parameter to determine the degree of multicollinearity between sets of observations (Liao and Valliant, 2012; García et al., 2015; Salmerón Gómez et al., 2016), and provides an indication of the effects of multicollinearity on the variance of the regression coefficients. A threshold of 10 is the most common value used as an indication of multicollinearity and a lower VIF is indicative of inconsequential multicollinearity (O’Brien, 2007). Multicollinearity was not considered problematic since the VIF values ranged from 1.38 to 6.11.

Model estimation was performed using the lavaan package (Rosseel, 2012) running in the R environment (R Development Core Team, 2011). The models were fitted with the maximum likelihood estimator, the most widely used in SEM. This estimator is unaffected by data transformation and has asymptotic properties, hence the minimum variance, unbiasedness, efficiency and consistency (Myung, 2003). The maximum likelihood estimator assumes multivariate normality but is appropriated and robust even if non-normality exists (Yuan and Bentler, 2005) and still produces centered estimates of the parameters (Finney and DiStefano, 2006) unless sk and ku are too severe (Kline, 2011), which was not the case for our data.

The model fit was assessed using the goodness-of-fit chi-square (χ²) statistic at a significance level α = 0.05. This model test statistic tests the null hypothesis that the population covariance matrix does not differ significantly from the model-implied covariance matrix. A failure to reject the null hypothesis means that the model is consistent with the sample data matrix, thus supporting the model assumptions. χ² is the most commonly used statistic for the model fit in SEM to appraise the overall model fit (Hu and Bentler, 1999) and it works particularly well for sample sizes between 75 to 200 cases (Kenny, 2015). It tends to be significant for larger sample sizes due to its sensitivity to sample size (Schumacker and Lomax, 2010). Additionally, the correlation residuals of the model fit were established to confirm the explanatory power of the model for the specific associations observed. As a rule of thumb, a correlation residual with an absolute value greater than 0.10 is a sign of a poor explanation of the corresponding sample correlation (Schumacker and Lomax, 2010). In addition to χ², several other approximate fit indexes were considered to provide further insights about the model fit such as the Steiger-Lind Root Mean Square Error of Approximation (RMSEA; Steiger, 1990), the Jöreskog-Sörbom Goodness of Fit Index (GFI; Jöreskog and Sörbom, 1982), the Bentler Comparative Fit Index (CFI, Bentler, 1990), the Tucker-Lewis Index (TLI; Tucker and Lewis, 1973), the Standardized Root Mean Square Residual (SRMR) and the ratio χ²/df. Taken together, these are among the most common approximate fix indexes that exist in the SEM literature (McDonald and Ho, 2002).

A multigroup analysis was also performed to determine how the considered disturbance factors specifically affect each riparian succession phase. This analysis was preceded by an evaluation of model measurement invariance, to judge the validity of the proposed comparisons, along with the computation of CFI variation (ΔCFI) to confirm the evaluation. During the multigroup analysis, the proposed models were individually fitted to each succession phase because the sample size of each group was considerably lower than the recommended sample size for these analyses. Consequently, a 95% confidence interval for the mean value of the group-specific model parameters was built in these cases based on larger sample sizes generated by the bootstrapping resampling method with 1000 bootstrap draws. All of the relationships between variables were analyzed using standardized coefficients for a better understanding of the direct effects between the variables (Grace and Bollen, 2005).

Figure 3 shows the proposed model construct for patch location, which was successfully fitted with a non-significant χ² test (p-value = 0.221). Accordingly, the population covariance matrix was not significantly different from the covariance matrix estimated by the model at a confidence level of 95%. The additional approximate fit indexes also supported this outcome and provided a perception of the goodness of the adjustment of the proposed model for the sampled data. The CFI, TLI, GFI and RMSEA indicated a very good adjustment, but the SRMR and χ²/df indicated that this result was a good adjustment (see Schumacker and Lomax, 2010; Kline, 2011; Marôco, 2014 for summaries on the fit classification indexes). Furthermore, all of the correlation residuals were less than 0.10, which suggested that the model could properly explain the corresponding sample correlations.

In this model, mrphd is manifested on MDi and grndh on GWDi. The latent variable of patch location is manifested on patch distance (THAD) and height (THAH) above thalweg, which indicates the patch positioning along the river stretches. The fitted model explains approximately 80% of the patch location variability, which implies that these variables contribute a considerable amount of information regarding patch location. The results show that the location of the riparian patches is primarily driven by grndh (standardized path coefficient = 0.930), whereas the mrphd seems of very little importance (standardized path coefficient = 0.061) to patch location. Regardless, grndh and mrphd are moderately correlated (-0.592). Consequently, riparian patches located at larger horizontal and vertical distances from the river thalweg appear to be subject to greater grndh disturbances. Furthermore, these results show that changes in grndh may lead to shifts in the location of succession phases, the eventual loss of the more vulnerable ones and the following extinction of its characteristic species.

The results for measurement invariance show that the patch location model indicates a weak invariance between study sites (H0: the model is invariant between study sites, not rejected for weak invariance with a p-value = 0.335). This finding means that the model construct is similar at all of the study sites and that for a given indicator, the factor loadings are significantly the same between the study sites. The ΔCFI (0.001) for the loadings is less than 0.01, indicating once more that invariance should not be rejected (Cheung and Rensvold, 2002). The model shows strong measurement invariance for the succession phases (the p-values for the weak and strong invariance tests are 0.443 and 0.167, respectively). The ΔCFI for the loadings and intercepts is correspondingly, 0.001 and 0.004, again smaller than the proposed threshold of 0.01 and therefore support the results of the invariance tests. A weak invariance allows the comparison of the relationships between the latent variables across groups while a strong invariance allows for the inter-group comparison of latent variable means and covariances (Chen, 2008).

The proposed model was successfully fitted in the multigroup analysis (Figure 4). The χ² test was not significant (p-value = 0.165) and the approximate fit indexes indicated the good adjustment of the model. The model explained more than 68% of the data variance, for group sample sizes of 74, 52, 72, and 22, respectively, for IP, PP, ES, and EF. As in the general location model, for every succession phase, both disturbances are still well correlated and the grndh is the main driver of patch location (the standardized path coefficients were greater than 0.81 for all succession phases). In contrast, the mrphd had a residual effect on the succession phases, except for ES (standardized path coefficient = 0.207), for which an increase in the mrphd determined greater distances to thalweg. This means that mrphd pushes away this succession phase to a more distant location, thus preventing vegetation encroachment. In detail, IP is characterized by recently disturbed patches where woody vegetation is starting to establish. The seedlings survive according to the recruitment box (Mahoney and Rood, 1998) and therefore, the link with the water table elevation is tight. If the recruitment survives, the patches evolve to PP after approximately 2 years. Individuals have now a settled root system but still vegetate inside the active channel, where the groundwater remains within reach almost all year round, accounting for the less pronounced link with the grndh. The ES patches live in a survival limbo between the grndh and mrphd. The indicator species of this succession phase (Salix spp.) are obligate phreatophytes well adapted to flow disturbances (González et al., 2010), whose existence is compelled by mrphd to the limit of the tolerated distance from the groundwater table. Consequently, these patches are considerably more likely to be negatively affected by groundwater level fluctuations and are therefore highly dependent on grndh. In turn, the typical facultative phreatophytes that characterize EF confer to this succession phase less dependency on this factor and a survival advantage over previous succession phases.

The 95% confidence intervals of the bootstrapped coefficients confirm the discrepancy in the effect between the grndh and mrphd on the location of the succession phase patches (Table 2). In fact, the magnitude of the effects of the two disturbances is completely distinct in the first three succession phases and always higher for grndh in all succession phases. Furthermore, one cannot rule out the possibility of that a mrphd would have no effect on the location of the succession phases, unlike grndh, which can be the only driving factor of the location of the succession phases as indicated by the proposed model. Succession phases have different heights and distances to thalweg, so it is possible to infer the location of the succession phases along the river stretch despite some degree of overlap. IP is the closest succession phase to the thalweg (the standardized values for mean distance and mean THAH is 3.220 and 3.406, respectively), followed by ES (4.669 and 4.731), PP (4.787 and 5.984) and EF in the outer parts of the river (5.499 and 6.824).

The patch shape model was successfully fitted with a non-significant χ² test (p-value = 0.078), and the additional approximate fit indexes indicated a good adjustment classification (Figure 5). GFI, CFI, TLI and SRMR indicated that this adjustment was very good while χ²/df and RMSEA classified it as acceptable. The latent variables of fluvial disturbance are manifested as GWDi and MDi while the latent variable patch shape is manifested as patch PERIMTR and the MPFD. Nevertheless, this model only explained approximately 13% of the patch shape variability and the disturbances accounted for have smaller effect on this feature (the sum of standardized patch coefficients is 0.562). The grndh continues to be the main driver of the two disturbance factors, despite that the effect of the mrphd (standardized path coefficient = -0.139) has now proportionally increased to approximately a third of the effect of the grndh (standardized path coefficient = -0.423).

The investigation of the measurement invariance revealed a weak and strong invariance between study sites (p-values of 0.268 and 0.844), respectively. Likewise, the ΔCFI for loadings (0.001) and intercepts (0.001) were smaller than the proposed threshold of 0.01, providing confirmation of both levels of measurement invariance. The model only shows weak measurement invariance (p-value = 0.383) for the succession phases supported by a ΔCFI for the loadings of 0.0004.

Notwithstanding, the patch shape model did not achieve a successful fit in the multigroup analysis. The χ² test was significant [χ²₍₈₎ = 20.791 and p-value = 0.008], meaning that the population covariance matrix was significantly different from the model-estimated covariance matrix. Also the approximate fit indexes indicated adjustment problems with contradictory classifications, specifically, an unacceptable fit by the TLI (0.757) and RMSEA (0.171), an acceptable fit by χ²/df (2.60), CFI (0.919) and SRMR (0.056), and a perfect fit by GFI (1.000). In addition, the explanatory power of the model was very low for all of the succession phases (0.066 for IP, 0.048 for PP, 0.034 for ES and 0.225 for EF). Taken together, these results demonstrated that the proposed model did not adequately explain the patch shape dynamics of the riparian succession phases in the Mediterranean rivers.