Drainage Network Patterns Determinism: A Comparison in Arid, Semi-Arid and Semi-Humid Area of Morocco Using Multifactorial Approach

Abstract

Drainage network patterns influence the hydrological response of the watersheds and must be taken into account in the management of the water resource. In this context, it is important to identify the factors that control the configuration of drainage networks in and beyond specific climatic conditions. Here, we study 318 basins spread over three sectors (arid, semi-arid, and semi-humid) of Morocco where seven drainage network patterns have been identified. From each basin, 14 parameters were extracted, describing the relief, geology, morphometry, drainage network, land cover, precipitation, and time of concentration (Tc). Principal component analysis (PCA) and discriminant analysis (DA) processing were performed on the entire database and on each sector separately. The results show that the drainage network pattern is a feature of the landscape that contributes significantly to the variance of the basins. They suggest that the distribution of network patterns is controlled by the relationship between the different parameters, mainly those related to the relief, more than by the variations of each parameter taken individually. The network discrimination rate is 63.8%, which improves when each sector is treated separately. Confusion in discrimination are similar across all sectors and can be explained by similar conditions (active tectonic, deformation, and uplift) or transitions from one network pattern to another, due to the landscape evolution of certain sectors. A contribution of climatic variables appears locally but was attributed to a statistical coincidence, these parameters presenting a distribution close to that of the relief and geology variables.

1. Introduction

Climate change tends to intensify the global hydrological cycle, in particular by increasing inequalities in the spatial and temporal distribution of precipitation and the frequency of extreme events [1]. The consequences are both longer periods of drought, modifying the climatic limits of arid regions, and floods, the severity of which seems to be increasing in many parts of the world [2]. It is now recognized that the drainage network pattern has a strong influence on the response of catchment areas to rainfall events [3]. Various regional or local constraints affect the development of the drainage network, leading to the formation of a wide variety of patterns with distinct characteristics and that resulted in a typology [4,5,6]. Previous works have suggested that the development of drainage network patterns is controlled by two main factors, namely slope and geological structuring [4,5,6]. These works lead to hypotheses and a classification of the drainage network patterns on the basis of findings, careful observations, and uni-factorial descriptions. However, the links between the physical properties of the environment and the type of drainage networks are neither ranked by importance nor strictly demonstrated, probably because this link is partial, as the determinism of the drainage network pattern can be more complex than a causal link based on a single parameter [7]. In particular, the relief parameters, climate aridity, and vegetation cover affect the extent of the dissection of the drainage basins. Precipitation, temperature, and water storage capacity are also relevant factors that control rock weathering and influence the formation and development of watercourses [8].

The combined approach of GIS (Geographic Information System) and multivariate statistical analysis seems suitable for the analysis of phenomena involving several physical characteristics and properties, as well as for the processing of a large amount of data [9], while searching for the sources of variability [10,11]. This is the approach we tested in this work, the aim of which is to understand the role of various environmental factors (relief, geology, basin morphometry, land cover, climate) in the formation and distribution of drainage network patterns. The study is based on classical statistical techniques such as principal component analysis (PCA) and discriminant analysis (DA) and on a set of 318 basins. In order to increase the geographical scope of the study, the basins studied are distributed over three sectors in Morocco, respectively arid, semi-arid and semi-humid.

2. Materials and Methods

2.1. Study Area

The External Rif, the Middle Atlas, and the Central Anti-Atlas were selected as their climate is classified as semi-humid, semi-arid, and arid, respectively (Figure 1a). These three sectors have also distinct and very complex geological and geomorphological conditions (Figure 1b). The External Rif, which is part of the Rif mountain range, is bordered by the Mediterranean Sea to the north and the Middle Atlas to the south. The western and central parts, exposed to both oceanic and Mediterranean disturbances, are very wet with up to 2000 mm of rainfall per year, which makes it the wettest region of the Maghreb. Further east, the climate is much warmer due to the eastern continental influence of hot and dry winds called ‘Chergui’. The External Rif consists mainly of high relief (Intrarif, Mesorif, and internal Prerif), and plains (External Prerif and Guercif plain) dominated by clay and shale formations, with progressive and rounded relief. Lastly, some limestone massifs are located in the northern part [12] (Figure 2a). The Middle Atlas refers to the vast mountainous area forming part of the alpine domain. It is bordered to the north by the South-Rifain corridor (Saiss plain), to the south by the High Atlas and the high Moulouya Valley, to the east by the middle Moulouya Valley and to the west by the Central Morocco. The western side of the Middle Atlas, in a situation of first ascendancy, receives on average 1000 to 1500 mm of rainfall per year, from oceanic weather disturbances. In contrast, the eastern valleys are much drier and impacted by desertification. This sector consists in a mountain range composed of volcanic plateaus and limestone rocks, namely the tabular Middle Atlas in the west, with an altitude of 1000 to 1500 m and the folded Middle Atlas in the northeast, with a higher altitude often exceeding 3000 m. Both are separated by a major fault called the North Middle Atlas Fault (NMAF) [13] (Figure 2b). The Anti-Atlas is the southernmost of the Atlas massifs and the driest mountain range in Morocco, receiving only 100 to 200 mm of rainfall per year. Water flows in some rare places, forming precious basins of clear water. It is an old eroded and desert low mountain range on the edge of the Sahara, with an average altitude of about 1500 m. Its center, the Central Anti-Atlas, includes several Precambrian massifs or “inliers” (The Siroua massif, The Bou Azzer-El Graara, The Zenaga inlier, and The Saghro massif) made up of Proterozoic terrains that are surmounted by a thick folded Palaeozoic sedimentary series [14]. The reliefs are incised with imposing gorges, and the mountain range is separated in two by the stitching of the Oued Drâa and the Major Anti-Atlas Accident (MAAA) (Figure 2c).

2.2. Selection of the Descriptive Parameters

The drainage basins and the drainage network were extracted from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) 30-metres resolution digital elevation models (DEM) using the flow accumulation method of Jenson and Domingue [17], and the analytical method of Tarboton et al. [18]. The drainage network has been identified with a threshold of 900 km², which was calibrated with respect to the heads of the real drainage channels obtained from 1/100,000 topographic maps. For each extracted basin, the geographical coordinates of the barycenter were determined, as well as 14 parameters that reflect the drainage network pattern, relief, geology, morphometry, land cover, and climate characteristics.

2.2.1. ‘Type’ of Network Parameter

Type refers the drainage network pattern in each basin. It was identified and assigned by visual interpretation according to several criteria: general shape of the network, orientation, tributary confluence angles, etc. (Figure 3) [4,5,6]. The following code was assigned to each pattern: dendritic = 0, pinnate = 1, parallel = 2, trellis =3, rectangular = 4, herringbone = 5 and barbed = 6. This order reflects an increasing structuring control on the drainage network [19].

2.2.2. Drainage Density

The parameter ‘Mean-Ddrain’ represents the average drainage density in each basin (Figure 4a). It is expressed as the ratio of total channel length to total basin area (km/km²) [20].

2.2.3. Morphometric Parameters

The morphometric characteristics of each basin in the three study areas were obtained from two widely used shape indices, which are easy to calculate and interpret. On the one hand, the Gravilus index “K_g” reflects the stretched shape of the watersheds, which gives an idea of the peak flow. K_g is high for an elongated shaped watershed and it is low for the almost circular shaped watershed.

$$K_g =\frac{P}{2\sqrt{\mathsf{\pi }}·A}$$

(1)

where “P” is the perimeter of the basin in kilometers.

The Horton index “K_h” of a watershed is the relationship between the surface area of the basin (km²) and the length of the main river in kilometers [21]. It is close to 1 for a watershed of a compact shape, and the lower the value below 1, the more the watershed is elongated.

$$K_h =\frac{A}{\mathrm{l}²}$$

(2)

where “l” is the total length of the main river and “A” is the surface of the basin in km².

2.2.4. Relief Parameters

“Slope” indicates the average slope of each basin that was calculated from the 30-m resolution DEM using the method of Burrough and Mcdonnell [22]. It consists of calculating for each cell the maximum rate of change in elevation values relative to neighboring cells (Figure 4b). “R1”, the relief ratio described by Melton [23], is defined as the ratio of the difference in elevation between the highest point and the outlet of the basin “r” and its perimeter “p”.

$$\mathrm{R}1=\frac{\left(ALTMax-ALTMin\right)}{P}=\frac{r}{p}$$

(3)

“R2”: Schumm’s dimensionless relief ratio “R2” [24] is defined as the ratio between “r” and the longest dimension of the basin parallel to the main drainage line “L” (Figure 4c).

$$\mathrm{R}2=\frac{r}{\mathrm{L}}$$

(4)

“L” is defined as the length of the equivalent rectangle:

$$\mathrm{L}=\frac{K_g\sqrt{A}}{1.12}\left[1+{\sqrt{1-\left(\frac{1.12}{K_g}\right)}}^2\right].$$

(5)

The roughness index “H” calculated as the product of the drainage density “Dd” (m/m²) and the basin relief “r” (m). It refers to the level of smoothness and roughness of the basin surface and its vulnerability to soil erosion [25].

2.2.5. Geological Parameters

The geological characteristics for each basin were extracted from different data sources for the three study areas, i.e., geological data of Africa provided by the United States Geological Survey (USGS, http://earthexplorer.usgs.gov/), Geological maps of Morocco at a scale of 1/1,000,000, 1/200,000, 1/100,000 and 1/50,000 [16,26,27,28,29,30,31], as well as several previous geological studies carried out on the three study areas. The geological characteristics of the basins are represented by three parameters: the ‘Rock-type’ related to the dominant rock type in each basin. For the Central Anti-Atlas, the rock types were grouped into four major groups, with the following coding: Quaternary sedimentary rock = 0, Palaeozoic sedimentary rock = 1, volcanic rock = 2 and granitic-meta-sedimentary rock = 3. However, for the Rif and the Middle Atlas areas, since their lithology is dominated by sedimentary rocks and particularly carbonate deposits, rock types were grouped as follows: Alluvial, Clay, and Silt = 0, Marl = 1, Conglomerate, Shale, and Flysch = 2, Sandstone, Dolomite, and Limestone = 3, volcanic rock = 4 and crystalline rocks = 5. This rock coding order roughly reflects a progression from weaker to harder rock types [32]. The second geological parameter reflects the lithological homogeneity or heterogeneity of the basin ‘Homo-Geo’. A basin covered over 75% of its surface by the same lithological formation was classified as homogeneous and assigned code 0; otherwise, it was classified as heterogeneous and assigned code 1 [32]. The third parameter, ‘Line-Dens’, reflects the density of the geological lineaments (km/km²) identified in each basin. Lineament identification was carried out using the spatial filtering method [33], which consists of applying a filter to the Landsat 8 OLI band 6 satellite image and the DEM hill-shade images of 30 m resolution for the three sectors. After merging of both sets of geological lineaments, duplicates were eliminated to avoid redundant information, and the geological maps were used for validation (Figure 5a).

2.2.6. Land Cover

The parameter ‘Landcover’ was extracted to determine the distribution of land cover in each basin, using the global land cover classification data layer produced annually from 2001 to 2013 Moderate-Resolution Imaging Spectro-radiometer (MODIS) Land Cover Type 1: International Geosphere-Biosphere Programme (IGBP) global vegetation classification scheme. The data are available at 500m resolution in standard MODIS grids (USGS), covering approximately 1200 km × 1200 km (≈10° × 10° at the equator) [34]. The following code was assigned to each land cover type: Forest, 0; Mixed forest, 1; Wooded savannah, 2; Grassland, 3; Savannah, 4; Shrub land, 5 and Bare soil with sparse vegetation, 6 (Figure 5b).

2.2.7. Time of Concentration

The parameter Tc (time of concentration) is used as an important parameter in the calculation of peak flows during floods [35]. We used the Natural Resources Conservation Service (NRCS) method where the time of concentration of a drainage network is estimated using the watershed lag, which is defined as a function of the hydraulic length of the watershed [35].

$$\mathrm{L}=\frac{l^{0.8}{\left(\mathrm{S}+1\right)}^{0.7}}{1900{\mathrm{Y}}^{0.5}}$$

(6)

where L is the propagation time of the basin (h), l is the length of the main channel (ft), Y is the average slope of the basin (percent), and S is the maximum potential retention (in). The maximum potential retention related to soil and land cover conditions represents the maximum potential retention of soil moisture and can be written as follows:

$$\mathrm{S}=(\frac{1000}{\mathrm{CN}}) -10$$

(7)

where CN is the curve number that represents the hydrological soil types. A low CN value reflects high soil infiltration rates and high CN value reflects a high soil runoff rates. CN was derived from the MODIS classification of land cover and hydrological soil groups (HSGs) at 250 m resolution under equable hydrological conditions [36]. Then, the concentration time is estimated as follows:

$$\mathrm{Tc}=\frac{\mathrm{L}}{0.6}$$

(8)

2.2.8. Precipitation Parameter

We used the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) to calculate the average annual precipitation over a 14-year period (2005–2019) for each sector (Figure 5c). CHIRPS is a 35+ year quasi-global rainfall dataset from 1981 to present, with a spanning range from 50°S to 50°N (and all longitudes). CHIRPS incorporates a 0.05° climate CHPclim satellite imagery, together with in situ station data to create a gridded rainfall time series for trend analysis and seasonal drought monitoring [37,38]. CHIRPS have negative biases (−2.01%) before 2000, while this underestimation was corrected after 2000. Temporal analysis and intensity distribution of the gauge-adjusted CHIRPS estimates are in line with the estimates of the Global Precipitation Climatology Centre (GPCC) in several regions of the world, such as the United States, Europe, and Africa [39].

2.3. Multivariate Analysis

A principal component analysis (PCA) was performed by diagonalization of the correlation matrix in order to identify and rank the different sources of variability within the extracted basins dataset. The principal components are linear combinations of the 14 parameters and thus behave as macro-parameters. They are orthogonal to each other and therefore represent independent sources of variability, i.e., independent associated processes. Taking into account the main PCs makes it possible to concentrate the information in a reduced number of factorial axis while losing a minimum of the information contained in the dataset, which constitutes a dimensional reduction of the data hyper-space [7,9]. To determine the impact of climate differences on the factors that control the variability in drainage network patterns, PCA was first applied to the entire dataset (318 basins) and then individually to each sector. Then, a discriminant analysis was performed to identify the combinations of parameters that best distinguish the different types of drainage network patterns. DA is a frequently used linear algebra method to discriminate multiple sets characterized by multiple parameters in a dataset. It consists in the search for a discriminant function that projects the points on a line passing through the centroids of the groups of scatter plots in an optimal direction so as to best discriminate the scatter plots.

3. Results

3.1. Drainage Network Distribution

The extracted 367 basins showed areas ranging from 6 to 619 km². Among them, 49 had an unidentified drainage network pattern, due to the presence of dams, drainage network deficiencies, or determination issues caused by the DEM resolution (

Table 1. Distribution of drainage network patterns in the three studied sectors.

	External Rif	Middle Atlas	Central Anti-Atlas
Dendritic	22	4	-
Pinnate	-	-	6
Parallel	13	36	20
Trellis	6	20	15
Rectangular	42	37	45
Herringbone	23	23	-
Barbed	3	3	-
Unidentified	17	13	19
Total	126	136	105

). These have been eliminated in the subsequent mathematical processing. In the External Rif and the Middle Atlas, six drainage network patterns types were identified, namely the dendritic, parallel, trellis, rectangular, herringbone, and barbed patterns (Figure 6a,b). The Central Anti-Atlas contains only four patterns: pinnate, parallel, trellis, and rectangular (Figure 6c).

3.2. Principal Component Analysis

3.2.1. Analysis of the Total Database (318 Observations × 14 Variables)

The two first principal components (PC1 and PC2) accounted for 30.6% and 16.4% of the total variance, respectively (

Table 2. Eigenvalues and percentage of explained variance by the first and second principal components for the total database and the three sectors.

	Total Database		External Rif		Middle Atlas		Central Anti-Atlas
	PC1	PC2	PC1	PC2	PC1	PC2	PC1	PC2
Eigenvalue	4.3	2.3	5.3	2	4.2	2.7	4.8	3.2
Variance (%)	30.6	16.4	37.6	14.3	29.7	19	34	22.5
Cum. %	30.6	47	37.6	51.9	29.7	48.7	34	56.5

). PC1 showed positive correlations with the drainage network pattern parameter “Type”, relief parameters (Slope, H, R1 and R2), geological parameter “Line-Dens” and morphometric parameter “K_h”. It was also negatively correlated with “K_g”, drainage density “Mean-Ddrain”, and “Tc” parameters (Figure 7a). PC2 showed positive correlations with the parameter “Type”, “Tc”, “H”, “Rock-type”, “Line-Dens”, and “K_g”. In addition, it presented negative correlations with the “Landcover” and “K_h” (Figure 7a).

3.2.2. Analysis of the Database for Each Sector

The first factorial plans (PC1-PC2) accounted for 51.9%, 48.7%, and 56.5% for the External Rif, the Middle Atlas, and The Central Anti-Atlas, respectively (Table 2). For the External Rif, PC1 (37.6% of the variance) was positively scored with the relief parameter “Slope, R1, R2 and H”, “Line-Dens”, “K_h” and “Type”, and negatively with “Mean-Ddrain”, “Tc” and “K_h” (Figure 7b). PC2 (14.3% of the variance) showed positive correlations with “K_g”, “H”, “Tc” and negative with “K_h” (Figure 7b).

In the Middle Atlas, PC1 explained almost 30% of the observed variance (Table 2) and was correlated with the same parameters as PC1 for the External Rif (Figure 7c). PC2 accounted for 19% of the variances, and was positively scored with “K_g” and negatively with “K_h”, just as for the External Rif. In addition, it was also positively correlated with “Slope” “Precipitation” and “Line-Dens”, and it was negatively correlated with “Mean-Ddrain”.

In the Central Anti-Atlas, PC1 accounted for 34% of the variance (Table 2). Just as for the other sectors, PC1 showed similar correlation with “Type”, “Slope”, “R1” and “R2”, “Line-Dens”, “Mean-Ddrain” and “Tc”, but it was also positively scored with “Rock-type” and negatively scored with “Landcover” and “Precipitation” (Figure 7d). PC2 explained 22.5%, of the variance, showing positive correlations with “R1” and “R2” and “K_h” and negative correlations with “K_g” and “Landcover”, i.e., PC2 showed more similarity with the PC2 obtained on all the data (Figure 7a).

3.3. Discriminant Analysis

3.3.1. Analysis of the Total Database

Six discriminant functions discriminated the drainage network patterns with a success rate of 63.8% (

Table 3. Confusion matrix of the discriminant analysis for the total database.

from\to	0	1	2	3	4	5	6	Total	% Correct
0	20	0	2	0	1	3	0	26	76.9%
1	0	6	0	0	0	0	0	6	100.0%
2	2	0	58	0	6	3	0	69	84.1%
3	0	0	7	8	11	15	0	41	19.5%
4	2	0	22	0	67	33	0	124	54.0%
5	1	0	3	0	4	38	0	46	82.6%
6	0	0	0	0	0	0	6	6	100.0%
Total	25	6	92	8	89	92	6	318	63.8%

). Five of the seven drainage network patterns had a classification rate higher than 75%, namely the dendritic, pinnate, parallel, herringbone, and barbed patterns. The trellis and rectangular patterns presented many confusions with classification rates of 19.5% and 54%, respectively (Table 3).

The first three discriminant functions accounted for 87.28% of the discrimination (

Table 4. Discriminant functions for the total database.

	F1	F2	F3	F4	F5	F6
Discrimination (%)	59.40	15.93	11.95	8.95	2.53	1.24
Cumulative %	59.40	75.33	87.28	96.23	98.76	100.0

). The first function F1 (59.4%), showed strong negative correlations with “Slope” and “Line-Dens” with values close to −0.71. It was positively scored with “Mean-Ddrain” (0.88) and secondarily with “Landcover” (0.3) (

Table 5. Correlations between variables and discriminant functions of the total database (higher values are in bold).

	F1	F2	F3
Slope	−0.71	0.05	−0.01
H	−0.10	0.17	0.20
R1	−0.24	−0.05	0.07
R2	−0.23	−0.02	0.12
Kg	0.19	0.42	0.18
Kh	−0.12	−0.27	0.01
Rock-type	−0.21	−0.22	0.46
Homo-Geo	−0.06	−0.01	−0.34
Line-Dens	−0.71	0.20	0.44
MEAN-Ddrain	0.88	0.18	0.20
Landcover	0.30	0.30	0.03
Tc	0.01	−0.18	0.55
Precipitation	−0.21	0.010	−0.66

). F2 (15.93%) was moderately correlated with “K_g” (0.42) and the “Landcover” (0.3). The third function F3 (11.95% of the discrimination), showed moderate positive correlations with the “Tc” (0.55), “Rock-type” (0.46), and “Line-Dens” (0.44) and negative correlations with “Precipitation” (−0.66) and “Homo-geo” (−0.34).

The observations in the F1-F2 plot (75.3% of the discrimination) are displayed in Figure 8. The axis F1 isolated dendritic, pinnate, and parallel patterns in the positive side of the plot, which had low scores of slope and lineaments density and high scores of drainage density and “Landcover”. By contrast, the trellis, rectangular, chevron, and barbed patterns were located in the negative side of F1, forming a compact cluster. The F2 axis made it possible to discriminate the pinnate pattern in the positive side, with higher “K_g” values and land cover types order, from the dendritic pattern in the negative side. The other patterns showed no significant discrimination along F2 (Figure 8) and F3 (Figure 9). The F3 axis discriminated the dendritic patterns in the negative scores followed by the herringbone, barbed, and pinnate pattern that were closer to zero. Parallel, trellis, and rectangular patterns were superimposed and poorly discriminated along F3.

3.3.2. Drainage Pattern Type Discrimination within Each Sector

The DA identified five discriminant functions that correctly distinguished 89% and 83.7% of the six types of drainage network patterns for the External Rif and Middle Atlas, respectively (

Table 6. Confusion matrix of the discriminant analysis for each sector. (Dendritic = 0, pinnate = 1, parallel = 2, trellis = 3, rectangular = 4, herringbone = 5 and barbed = 6).

	From/to	0	1	2	3	4	5	6	Total	% Correct
External Rif	0	20	_	0	0	2	0	0	22	90.9%
	2	0	_	13	0	0	0	0	13	100%
	3	0	_	0	6	0	0	0	6	100%
	4	4	_	0	0	35	3	0	42	83.3%
	5	1	_	0	0	2	20	0	23	86.9%
	6	0	_	0	0	0	0	3	3	100%
	Total	25	_	13	6	39	23	3	109	89%
Middle Atlas	0	4	_	0	0	0	0	0	4	100%
	2	0	_	32	1	3	0	0	36	88.9%
	3	0	_	0	19	1	0	0	20	95%
	4	0	_	3	6	26	2	0	37	70.3%
	5	0	_	0	1	3	19	0	23	82.6%
	6	0	_	0	0	0	0	3	3	100%
	Total	4	_	35	27	33	21	3	123	83.7%
CAA	1	_	6	0	0	0	_	_	6	100%
	2	_	0	20	0	0	_	_	20	100%
	3	_	0	3	12	0	_	_	15	80%
	4	_	0	9	1	35	_	_	45	77.8%
	Total		6	32	13	35	_	_	86	84.9%

and

Table 7. Discriminant functions and their percentage of classifications for each sector.

		F1	F2	F3	F4	F5
External Rif	Discrimination (%)	58.90	17.84	14.12	5.64	3.50
	Cum. %	58.90	76.74	90.86	96.50	100
Middle Atlas	Discrimination (%)	66.49	13.05	10.39	5.83	4.24
	Cumulative %	66.49	79.54	89.93	95.76	100
Central Anti Atlas	Discrimination (%)	66.54	25.73	7.73	_	_
	Cumulative %	66.54	92.27	100	_	_

). In the Central Anti-Atlas, three functions discriminated among the four drainage network patterns with a success rate of 84.88% (Table 6 and Table 7).

In the External Rif sector, the drainage network patterns have been discriminated with a high success rate, but some residual confusion is observed for the rectangular, herringbone, and dendritic patterns (Table 6). F1 (58.90% of the discrimination) showed high negative correlations with “Slope” (−0.73), “Line-Dens” (−0.58), and the “Precipitation” (−0.46) and positive correlations with “Mean-Drain” (0.88) (

Table 8. Correlations between variables and discriminant functions of the three study areas (higher values are in bold).

Parameters	External Rif			Middle Atlas			Central Anti-Atlas
	F1	F2	F3	F1	F2	F3	F1	F2	F3
Slope	−0.73	−0.15	0.34	0.72	0.13	0.32	0.78	−0.05	0.10
H	−0.17	−0.39	0.34	0.08	0.1	0.38	0.22	0.64	0.15
R1	−0.21	−0.1	0.49	0.19	0.16	0.28	0.36	0.46	−0.04
R2	−0.21	−0.07	0.49	0.18	0.15	0.29	0.35	0.45	−0.04
Kg	0.19	−0.58	−0.07	−0.04	0.08	0.02	−0.25	−0.21	0.19
Kh	−0.1	0.54	0.17	0.03	−0.12	0.17	0.22	0.23	−0.18
Rock-type	−0.01	0.05	0.05	−0.06	−0.3	−0.1	0.78	−0.22	0.38
Homo-Geo	−0.07	−0.07	−0.39	0.13	0.10	−0.09	0.02	0.36	−0.01
Line-Dens	−0.58	−0.02	0.66	0.88	0.09	−0.23	0.65	−0.16	−0.09
Mean-Ddrain	0.88	−0.13	−0.05	−0.80	0.35	0.07	−0.67	0.42	0.31
Landcover	−0.06	0.53	−0.16	−0.11	0.58	−0.08	−0.51	0.35	−0.44
Tc	0.35	−0.16	−0.31	−0.22	−0.37	−0.02	0.55	−0.15	0.32
Precipitation	−0.46	−0.1	−0.25	0.75	0.35	−0.12	−0.8	0.1	0.10

). The F1 axis isolated the trellis, herringbone and barbed patterns in the negative side (Figure 10a), suggesting basins with a steep slope, high lineament density, and higher precipitation rates, by contrast with dendritic and parallel pattern with centroids located in the positive side. F2 (17.84%) showed positive correlations with “K_h” (0.54) and “Landcover” (0.53) and negative correlations with “K_g” (−0.58) (Table 8). The function F2 separated dendritic and barbed patterns in the positive side, i.e., wide-shaped basins with low vegetation cover, while trellis and parallel patterns are isolated in the negative side, suggesting elongated basins with more developed vegetation cover (Figure 10a). F3 (14.12%) presented positive correlation values for “Line-Dens” (0.66), “R1” and “R2” (0.49), “H” and “Slope” (0.34) and negative values for “Homo-Geo” (−0.39) and Tc (−0.31). The F3 axis (Figure 10b) made it possible to discriminate clearly the barbed pattern, whose centroid showed very high coordinates on the positive side, suggesting that this pattern occurs in conditions of rugged and high relief with high lineament densities. The centroids of trellis and dendritic patterns presented negative values on F3, indicating some lithological heterogeneity and/or high Tc. The centroids of herringbone, rectangular, and parallel patterns are close to zero.

In the Middle Atlas sector, all the drainage network patterns presented high discrimination rates, with a 100% discrimination score for the dendritic and barbed patterns. The lowest rate (70.3%) was for the rectangular pattern with 11 misclassified basins among 37 (Table 6). F1 (66.49%) showed positive correlations with “Line-Dens” (0.88), “Slope” (0.71) and the “Precipitation” (0.75) and negative correlations with “Mean-Ddrain” (−0.80) (Table 8). Just as for the previous sector, the F1 axis separated trellis, herringbone and barbed patterns, from the dendritic and parallel patterns, suggesting similar conditions for the development of these patterns (Figure 10c). F2 (13.05%) presented positive correlations values for “Landcover” (0.58) and “Precipitation” (0.35) and negative correlations for “Tc” (−0.37). The centroids of dendritic and barbed patterns were located on the negative side of F2 axis. F3 (10.39%) showed positive correlations with “Slope” (0.32), and “H” (0.38). It separated herringbone and barbed patterns with positive values from other patterns with centroids coordinates near-zero or negative (Figure 10d).

Finally, in the Central Anti-Atlas, pinnate and parallel patterns were discriminated at 100%, and the trellis pattern displayed three confusions with a discrimination rate of 80%. Rectangular patterns showed the weakest discriminations rate with 10 confusions among 45 basins (Table 6). F1 (66.54%) showed high positive correlations with “Slope” and “Rock-type” (0.78), “Line-Dens” (0.65), and “Tc” (0.55) and negative correlations with “Precipitation” (−0.8), “Mean-Ddrain” (−0.67), and the “Landcover” (−0.51). From these conditions, the F1 axis allowed isolating rectangular and trellis patterns that formed an isolated cluster with positive coordinates (Figure 10e), suggesting steep slope conditions, with predominantly granitic–metamorphic rock type, with high lineaments densities, low drainage densities, a relatively developed vegetation, and a high Tc. Pinnate patterns were well discriminated on the far negative side. Both parallel and pinnate patterns showed negative coordinates, indicating that they developed in conditions opposite to those described for rectangular and trellis patterns. F2 (25.73%) showed positive correlation with “H” (0.64), “R1” (0.46), and “R2” (0.45), followed by “Mean-Ddrain” (0.42). A parallel pattern appeared on the positive side of F2, suggesting basins that occur in fairly high relief with a high drainage density. F2 placed the pinnate pattern in the negative side, suggesting low-relief conditions. Rectangular and trellis patterns presented centroids with coordinates in the negative side of this axis, but they were weakly discriminated. F3 (7.73%) presented a moderate positive value for “Rock-type” (0.38) and a negative value for “Landcover” (0.44). It discriminated rectangular and trellis patterns that were almost superimposed on the F1 and F2 axis (Figure 10f).

4. Discussion

4.1. Multifactorial Influence in Determining the Type of Drainage Network

The first factorial axis (PC1) is correlated with the same parameters, whether we consider all the data or each sector separately. This shows the importance of PC1 as a macro parameter beyond regional heterogeneity. The associated eigenvalue is 4.2 to 5.3 (Table 2), depending on the processing, which means that it represents as much variability as four to five variables, which is considerable. We also note a strong contribution of the drainage network pattern to PC1 (Figure 7). The same is true for the results of the discriminant analysis. The first function F1 that displays a discrimination rate of 63.8% for 318 basins and seven drainage network patterns is very satisfactory. F1 shows high correlations with the same parameters on each of the three sectors (Table 8) and on the entire database (Table 5). It mainly depends on the slope, geological structure, drainage density, and land cover. Therefore, the drainage network pattern is a landscape characteristic that strongly contributes to the variability of the studied basins, in correlation with topographical, geological, and climatic factors. The previous works [4,5,19] demonstrated, from single-factor analyses and numerous observations, that the slope and geological structure were the major factors controlling the distribution of drainage network patterns. The correlation between these parameters in this study highlights the multifactorial aspect of the process of drainage network patterns determination, thus suggesting that it is the relationship between the different topographic, geological, and climatic parameters that controls the distribution of each drainage network pattern and not the variability of each parameter taken individually. Considering each sector individually, the discrimination rate increase slightly, but the order of discrimination of the network patterns remains similar (Figure 10a,c,e). Therefore, the major multifactorial process that controls the formation of drainage network patterns is general and depends little on local conditions of each sector.

4.2. Details of the Sources of Variability

4.2.1. Drainage Pattern, Relief, and Geological Parameters

The positive correlation between “type” and the four relief parameters (slope, R1, R2, and H) (Figure 7) on PC1 suggests that the relief conditions largely control the distribution of the drainage network patterns. This can be explained by the two families of drainage pattern developed in our study area [40]. On the one hand, the herringbone and barbed patterns are high relief patterns, generally developing in high areas with a strong drop and lateral slopes higher than the regional slope, all combined with strong tectonic activity and intense deformations. These conditions are well represented in the sectors of the External Rif and the Middle Atlas [41,42]. On the other hand, the dendritic, pinnate, parallel, trellis, and rectangular patterns are patterns of low relief, that is to say of a weakly eroded zone associated with a weak tectonic activity, as is the case in the Central Anti Atlas [14]. The correlation between PC1 and the geological parameters “Line-Dens” and “Rock-type” is linked to the local structural control, highlighted by the presence of drainage network patterns (rectangular, trellis, herringbone, and barbed) whose development is mainly controlled by geological structuring [4,5,7,19,23,43].

Similarly, F1 separates two groups of drainage network models. The first group consists of patterns controlled by rugged relief and strong structural control (barbed, herringbone, rectangular, and trellis) (Figure 8d–g). These patterns develop in structural areas with a high density of geological lineaments and deformations due to intense uplift tectonic activity [7,42,43,44]. The uprising concerns the Intrarif, the Mesorif, and the internal Prerif in the External Rif sector, and the folded and tabular Middle Atlas, central Morocco, the Paleozoic bedrock in the Middle Atlas sector, and the Saghro, Bou Azzer, Siroua, and Agdez massifs in the central Anti-Atlas sector (Figure 11). The second group includes parallel, pinnate, and dendritic patterns. These are located in regions of low topographic amplitude and weak geological structure, namely the plain of Guercif and the external Prerif for the External Rif sector, the Saiss and Missour plains in the Middle Atlas sector, and the Feijas escarpment and its alluvial plains in the central Anti-Atlas sector (Figure 12).

As mentioned above, regions where barbed, herringbone, trellis, and rectangular patterns predominate exhibit rugged relief but also several areas of weakness and associated tectonic structures that limit the connection of lower-order flows, thus preventing the further development of immediate higher-order flows. This leads to a decrease in both the total length and the number of streams, and it ultimately results in low drainage density [10,19,43], which explains the negative correlation of these patterns with the drainage density parameter. In the low relief areas, with low lineament density, the drainage network is controlled by slope rather than influenced by structural control, and pinnate and parallel patterns prevail. Coupled with the prevalence of sedimentary rocks, it leads to a high drainage density [7,19,20,23,45].

4.2.2. Drainage Network and Basin Morphometry

The correlation of the morphometric parameters K_g and K_h with PC1 and PC2 (Figure 7a) obtained on all the data reflects the importance of morphometry in the variability of basins and its relationship with the drainage network pattern. However, this contribution is not general over the entire study area, since it does not appear in the PC1 of the Central Anti-Atlas sector (Figure 7d). The External Rif and the Middle Atlas sectors are geologically recent structural domains compared to the Anti-Atlas, which belongs to a more stable African craton. In these two sectors, the basins present a fairly young stage of erosion due to active alpine tectonics and deformations, in particular uplifts, which characterize the evolution of the landscape and the distribution of drainage networks [44,46,47]. The high correlation scores observed with the morphometric, relief, and lineament density parameters corroborate this interpretation (Figure 7b,c). More humid climates are often confined to mountainous regions [23], that is to say often on resistant rocks affected by significant tectonic activity, and they are likely to show a correlation between morphometric and precipitation parameters during statistical analyses, as observed in the Middle Atlas sector (Figure 7c).

4.2.3. Drainage Networks and Hydrological Response

Taking into account the whole database, the time of concentration “Tc” is correlated with the first two principal components, suggesting that the drainage network pattern influences the hydrological response of the basins. The correlation is negative on the first factorial axis (PC1) (Figure 7a), indicating that the drainage network patterns associated with high landforms and high structural control have lower Tc, i.e., a faster hydrological response. However, this correlation being distributed on two factorial axes (PC1 and PC2) orthogonal in the hyperspace of the data, it must be attributed to distinct and independent processes. Both the External Rif and the Middle Atlas are part of the alpine domain characterized by the highest altitudes in Morocco and active uplift tectonics. These conditions generally coincide with strong river power [48]. In the central Anti-Atlas, basins with rectangular and trellis patterns have higher Tc than basins with parallel or pinnate patterns, which are associated with bare land or very low vegetation cover, steep slopes, narrow, and elongated shapes—that is to say, a set of conditions favorable to runoff. In other words, when processing all the data, PC1 reflects the hydrological response of the basins according to the conditions of the semi-humid External Rif and the semi-arid Middle Atlas, while PC2 reflects the arid conditions of the Central Anti-Atlas.

The results of the discriminant analysis (F1, Table 8) suggest that the time of concentration “Tc” contributes to the discrimination of the drainage network patterns only in the central Anti-Atlas sector. This result is in agreement with several studies that mention the influence of the drainage network patterns on the flooding risks sensitivity of arid zones in Morocco or in other regions of the world [35,49,50].

4.2.4. Climate Variability and Types of Drainage Networks

The “Land cover” and “precipitation” parameters contribute to the first factorial axis (PC1) over the central Anti-Atlas sector (Figure 7d), which at first glance seems surprising due to its extreme aridity. “Landcover” also contributes significantly to the F1 on this sector. This finding does not necessarily mean that the distribution of the drainage network patterns is directly influenced by climate parameters. Indeed, the land cover map of the central Anti-Atlas (Figure 5b) shows a strong similarity with the distribution of the drainage networks patterns. To the west, north, and north-west, the basins covered with open shrub vegetation of the arid type present mainly rectangular and trellis patterns. On the other hand, to the east, south, and south-east of the sector, parallel and pinnate patterns prevail, in basins with bare land or very sparse vegetation. On the one hand, the explanation is in the exposure of the drainage network patterns that develop at the massifs and high altitudes, such as the rectangular and trellis patterns, to Atlantic maritime influences, thus receiving more precipitations. These relief and climate conditions, which are associated with the fractured terrain containing deep valleys with a more humid climate and water springs usually located close to the tectonic accidents, lead to the development of the vegetation (Figure 13a). On the other hand, moving to the south-east of these areas, the altitude gradually decreases, together with a climate shift from semi-arid to extremely arid due to the hot and dry Saharan winds mainly in the plains and stony plateaus of Hamadas (Figure 13b) with sparse arid vegetation or bare land [49]. Therefore, the statistical correlation observed between the climatic parameters and the drainage network patterns, which also appears as a contribution to the discriminant function F2 (Table 5), corresponds to a control by the geological and relief parameters, these having the same distribution as the climatic parameters, a phenomenon already described by several authors [19,23,51].

4.3. Discrimination and Confusion in the Discrimination of Network Patterns

It is difficult to distinguish a process associated with discriminant functions that account only for a small percentage of discrimination, as in the case of F3. However, it appears that on each sector, F3 allows the discrimination of a drainage network pattern. On the one hand, the patterns of dendritic and pinnate networks are discriminated on the basis of the geological parameters “Rock-type” and “Line-Dens” because they develop on sedimentary and/or quaternary surface rock types in the absence of prominent tectonic structures (Figure 12). On the other hand, the herringbone and barbed patterns are associated to the “Tc”, “Precipitation”, and “Homo-Geo” parameters, as they occur in high relief of the External Rif and Middle Atlas sectors are exposed to humid Atlantic influences [52] with a faster basin hydrological response. F3 shows a negative correlation with the “Homo-Geo” parameter (Table 5). This reflects a lithological heterogeneity and thus confirms the conditions of a complex geological structure, in which these two high relief patterns evolve [40,42,47].

The main confusions in discrimination concern rectangular and trellis patterns, which are confused with each other, or with the parallel pattern (Table 3). This type of confusion was expected, because both rectangular and trellis patterns are generated by geological structuring [4,5,23], which are here developed on the same structural areas in the three sectors and connected to each other (Figure 6). Therefore, they are subject to the same conditions of relief and geological structure. In addition, the rectangular, trellis, and parallel patterns are part of the basic patterns likely to undergo the transition phenomenon, i.e., the passage from one pattern to another in space or time. It is generally due to the tectonic activity of uplift or subsidence, and in our case more specifically of uplift, as mentioned above [5,7].

Dendritic and pinnate patterns develop under similar conditions, and it is recognized that the dendritic type can evolve into a pinnate pattern when the basin consists of easily erodible material [4,5]. F2 discriminates these two patterns from “K_g” and “Landcover” parameters. The dendritic pattern concerns less elongated basins, which is explained by the low relief and poor conditions of structural control, leading to drainage networks with irregular ramifications in all directions, and tributaries joining the main course under all angles (Figure 12a–c), giving an often wide apple-tree-shaped basin [4]. The pinnate-pattern basins are only observed in the south-east of the Central Anti-Atlas sector, i.e., in conditions of extreme aridity and therefore with a high value of the Landcover parameter (Figure 12d). In this context, many almost parallel and closely spaced tributaries originate in the escarpment of semi-tabular monocline structures called Bani [45]. They join the main streams at an acute angle generating a feather-shaped drainage network (Figure 6c). The rate of discrimination increases slightly when considering each sector separately (Table 6), which reveals that drainage network patterns are better discriminated in the same environment, where they share the same topographical, geological, and climatic context. However, confusions remain between the same patterns (Table 3 and Table 6). Therefore, such confusions are general and do not depend on specific regional conditions.

5. Conclusions

This work is based on a dataset of 318 basins and 14 parameters related to the drainage network pattern, relief, morphometry, geology, land cover, climate, and hydrological response. It aimed to understand, by conventional statistical methods, the factors that control the development and distribution of drainage network patters in three sectors of Morocco in a climatic range from arid to semi-humid. Seven drainage network patterns were identified, namely Herringbone, barbed, rectangular, trellis, parallel, pinnate, and dendritic. The results of the principal component analysis confirm that the configuration of the drainage network explains an important part of the variance of the basins, while being correlated with the relief (slope, roughness, and texture) the geological conditions (type of rock and geological structure), the drainage density, and the hydrological response, over the entire study area or over each sector analyzed individually. This suggests that it is the relationship between the different parameters that controls the distribution of each drainage network pattern and not the variability of each parameter taken individually. A contribution of climatic variables appears locally but indirectly through the interrelation with the relief and geology variables over the study area. The rate of discrimination among drainage network patterns is in the order of 64%, reaching over 80%, when each climatic sector is treated individually. Confusions between patterns can be explained by similar formation conditions or by transitions from one pattern to another. These confusions are similar whether we consider the entire database or each sector. Therefore, they are general and do not depend on specific regional conditions. The influence of climate on the formation and distribution of drainage networks should be studied over a wider climatic range. This work was based on multivariate statistical methods using linear algebra; i.e., it is assumed that the relationships between the different descriptors are linear, which is not necessarily the case for a natural environment. Non-linear methods, such as machine learning, which can provide a gain in the discrimination of network types, should be tested in order to check whether the hypothesis of linearity corresponds to a real constraint. Finally, as the climatic range taken into consideration in this study is limited, such a study could be attempted on environments that are much more contrasted in terms of climate and plant cover, with the aim of establishing a work with a wide geographical scope.