Tarik Bouramtane

Suzanne Yameogo

Meryem Touzani

Abdessamad Tiouiouine

M 'hamed

El Janati

Jamila Ouardi

Ilias Kacimi

Vincent Valles

Laurent Barbiero

Statistical approach of factors controlling drainage network patterns in arid areas. Application to the Eastern Anti Atlas (Morocco)

Keywords: Drainage network patterns, Arid area, Principal Component Analysis, Linear discriminant analysis, Anti-Atlas, Morocco. 1. Introduction

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

geology and drainage pattern considerations are frequently used to interpret the geological context without previous mapping [START_REF] Drury | Image interpretation in Geology, 3rd editio[END_REF]. The differences in drainage patterns have led geomorphologists to develop classification systems for channel networks and to extract clues and information that reinforce hypothesis on some major geological events. [START_REF] Johnson | Streams and Their Significance[END_REF] conducted the first work on the meaning of drainage patterns, assuming that rivers follow the slope first and then adjust to structure. [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF] observed that surface slope and rock resistance determine drainage patterns. This author stated that slope induces the formation of patterns such as dendritic, parallel, radial and distributary, while geological structures (Joints, Faults and folds) produces straight, angular, trellis and annular arrangements. [START_REF] Parvis | Drainage pattern significance in airphoto identification of soils and bedrocks[END_REF] expanded this classification system by introducing fifteen modified patterns. [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF] conducted a more comprehensive study, adding new basic patterns and giving them a broader definition. He also suggested that transition patterns could result from changes from one scheme to another over time.

The eastern Anti-Atlas of Morocco (EAA) is characterized by a wide variety of prominent and distinct geological structures, ranging in age from Precambrian to Quaternary. The EAA have suffered several geological events and constraints of various orientations, such as the Hercynian collision between Laurentia and Gondwana responsible for the reactivation of several inherited faults, a rise of the Precambrian basement and a folding of the Paleozoic cover [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF]. In this context, deciphering the role of the parameters controlling the development of the various types of drainage network necessarily involves a statistical approach capable of extracting information contained in large databases that combines numerous observations and their characteristics. The objective of this work, carried out in the little known but complex area of the EAA, is to apply and combine various mathematical approaches to study this relationship between the different drainage networks patterns, topography and geology. A large number of parameters is used for the characterization of the area. An approach on individually taken or crossreferenced parameters 2 to 2 is quickly presented, but the complexity of the dataset requires a multivariate statistical approach. The independent sources of variability are studied by principal component analysis (PCA), which allows the information to be concentrated in a smaller number of macroparameters, i.e. the Principal Components. On this basis, hierarchical agglomerative clustering (HAC) allows a comparison of different types of drainage networks. Finally (and especially) a multiparameter linear discriminant analysis (LDA) makes it possible to determine whether it is possible to predict the type of drainage network based on the geological and topographical characteristics of the area, i.e. to verify whether there is a relationship between the environmental characteristics and the type of drainage network.

Material and methods

Study area:

The studied basins are located in the Eastern Anti-Atlas, in southeastern Morocco (Fig. 1). The Anti-Atlas forms a vast anticlinal extending from the Atlantic Ocean South-Westward to Hamada of Guir North-Eastward. The EAA represents sectors of Neoproterozoic paleo continents preserved from the Pan African to Alpine orogenies as recorded in their marginal belts. Intense foldings and overthrustings mark the limits between the Neoproterozoic belts and the surrounding cover (Fig. 1). The EAA region consists of six identifiable sectors (Table 1):

(1)

The Bani sector (BS) is a folded Paleozoic cover located southwest of the study area. The term "bani" refers to monoclinal semi-tubular structures from Cambrian to Devonian [START_REF] Choubert | Sur le Géorgien de l'Anti-Atlas[END_REF][START_REF] Colmenar | A synopsis of Late Ordovician brachiopod diversity in the Anti-Atlas, Morocco[END_REF][START_REF] Ouanaimi | Unraveling a distal segment of the West African Craton Paleozoic margin: Stratigraphy of the Mougueur inlier of the eastern High Atlas, Morocco[END_REF]. The Ordovician covers a vast area, with two ridges of hard sandstone, namely the first and second bani, separated by two levels of clay, called Feija and Ktaoua shales, respectively.

(2)

The sector of Ougnat-Ouzina axis (OOA) also is referred as the anticlinorium of Taouz in the literature. It is a wide NNW-SSE trend shear zone containing Cambrian anticlinal folded structures (Jbel Taklimt, Jbel Renneg and Jbel Tijekht). Folded Ordovician structural units are also described in the northeastern and southeastern part of this sector (Bou Mays and Amelane-Mech Irdane, the anticlinal units of Shayb Arras and Jbel Tazout), and are separated from each other by synclinal folds of Devonian and Carboniferous age (Ottara and Amessoui synclinals) [START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF].

(3)

The Maider Basin sector (MB), around the Lower Carboniferous Fesou synclinal in the center of the EAA is weakly deformed compared to other sectors. Its western boundary is a monoclinal folds zone, whereas in the south it is delimited by anticline folds trending E-W to ENE-WSW [START_REF] Kaufmann | Middle Devonian reef and mud mounds on a carbonate ramp: Mader Basin (eastern Anti-Atlas, Morocco)[END_REF][START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF].

(4)

The Kem-Kem Domain (KKD) includes a large stony tabular area of Cretaceous deposits, which highlights the partial erosion of sedimentary cover formed by detrital deposits. These formations remain intact, without significant developed structures [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF].

(5)

The South-Tafilalt Basin sector (STB) is characterized by the broad Lower Carboniferous Merzouga synclinal with a WNW-ESE trend axis, and parallel to the middle Devonian Erfoud anticlinal located northeast the sector.

(6)

The sector of the Saghro-Ougnat axis (SOA) comprises two large inliers, the Saghro and Ougnat Precambrian massifs, which constitute the highest points in the study area, at over 3200 m altitude. This sector is considered as a Pan-African arc [START_REF] Walsh | Neoproterozoic tectonic evolution of the Jebel Saghro and Bou Azzer-El Graara inliers, eastern and central Anti-Atlas, Morocco[END_REF], separating two significant geological domains, namely the High Atlas in the north and the Sahara of Morocco in the south. The inliers overlie deep faults reactivated on a scale of several kilometres, namely Assif n'Oussif Fault, Boulghzazil-Tinifift Fault, Akrouz Fault and Tizi n'Ressas Fault [START_REF] Soulaimani | Late Ediacaran-Cambrian structures and their reactivation during the Variscan and Alpine cycles in the Anti-Atlas (Morocco)[END_REF]. Several volcanic units occur and consist of andesitic to rhyo-dacitic volcanic and volcaniclastic rocks including mainly ignimbrites and tuffs. BS, MB and STB are lowland areas alternating with Appalachian landforms and mountainous areas in the north. MB and STB are more precisely plains of a synclinal nature. These lowland areas contain basin structures such as the Tazzarine basin in SB, the large fezzou syncline in MB and the Ziz plain in STB with quaternary deposits. Fig. 1. Geological context of the studied area and location of the 6 geological sectors Modified from [START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF]. Precambrian basement that outcrops in two large inliers, extremely fractured by important tectonic structures. This analysis was carried out on 230 basins extracted from 5 major endorheic basins, namely from the northeast to the southwest, the Ziz, Rhéris, Kem-kem, Maider and Draa basins. These endorheic basins and the procedure for the basins and channel network extraction from the ASTER-Global Digital Elevation Model are presented in Supplementary Material S1 and S2. The climate is arid in the north in both Saghro and Ougnat mountains, with annual rainfall of about 120 mm, and becomes extremely arid in the south, with only 50 mm per year [START_REF] Hilali | Hydrogéologie et ressources en eau du Tafilalet et ses régions limitrophes (sud-est du Maroc): connaissance, prospection, caractérisation, exploitation et gestion des ressources en eau[END_REF][START_REF] Kassou | Hydrogéologie, hydrogéochimie et qualité des eaux des nappes aquifères du bassin de l'Oued Reg, région d'Alnif[END_REF]. 2.2. Major parameters assigned to the basins: For each of the extracted basin, the geographical coordinates of the barycenter were determined, as well as 16 parameters that reflect their topographic, geological and morphometric characteristics. These parameters were obtained through various remote sensing data (Landsat 8 OLI, Sentinel-1A and Sentinel-2A), a Digital Elevation Model (30 m resolution ASTER-GDEM), 5 geological maps of 1/200000 scale that cover the entire study area [START_REF] Hinderleyer | Carte Geologique du Maroc au 1/200 000, feuille Jbel Saghro-Dades[END_REF][START_REF] Clariond | Carte Geologique du Maroc au 1/200 000, feuille BOU Haiara-Zegdou[END_REF][START_REF] Destombes | Carte Geologique du Maroc au 1/200 000, feuille TafilalteTaouz[END_REF][START_REF] Du Dresnay | Carte Geologique du Maroc au 1/200 000, feuille Todrgha-Ma[END_REF][START_REF] Choubert | Carte Geologique du Maroc au 1/200 000, feuille Zagora-Coude du Dra-Hamada du Dra[END_REF], and geological data of eastern Anti-Atlas from the United State Geological Survey (USGS, http://earthexplorer.usgs.gov/). The 16 parameters are detailed below.

Drainage network patterns parameter

The first is the type of drainage pattern identified in the extracted basins, with a code assigned to each pattern, as follows: 0 = dendritic, 1 = parallel, 2 = trellis, 3 = angular. The first step was to verify, for each basin, the adequacy between the extracted drainage network and the geological (scale 1/200000) and topographic (scale 1/100000) maps. Then the drainage pattern parameter was assigned on several criteria (general shape of the network, orientation, connection angles of tributaries, etc.) for recognition of drainage network patterns [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF][START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. In order to investigate the change in the drainage network orientation in the study area, the drainage network has been classified according to the associated rock type and structural regime. Five drainage network groups were identified, developed on the Precambrian bedrock, Paleozoic outcrops, Lower Mesozoic formations, Lower Neogene and those formed in Quaternary deposits.

Morphometric parameters

Morphometric characteristics for each drainage basin were calculated, namely the total area, the perimeter, as well as two shape indices known as the Gravilus compacity index (KG) that reflects the stretched form of the catchments, and the Horton shape index (KH), wish is indicative of the flood regime of the stream [START_REF] Horton | Drainage-basin characteristics[END_REF]. KG is defined as the ratio between the perimeter (P) of the basin to the perimeter of the circle having the same area (A):

𝐾 𝐺 = 𝑃 2√π.𝐴

(1) KH expresses the ratio of the basin area (A) to the length of the mainstream (L):

𝐾 𝐻 = 𝐴 L²

(2) 2.2.3. Topographic parameters In order to identify which feature of the relief (Fig. 2a) best explains the hydrological behavior of the basins and the structure of the drainage network, 6 parameters were calculated. They are: 1 -MAX-Alt is the maximum elevation above the basin outlet, used to express the basin relief [START_REF] Sangireddy | Controls of climate, topography, vegetation, and lithology on drainage density extracted from high resolution topography data[END_REF]. 2 -R1 is the relief ratio described by [START_REF] Melton | An analysis of the relations among elements of climate, surface properties, and geomorphology[END_REF], and defined as the ratio of basin 'r' wish is the relief difference in altitude between the culminating point and the outlet, and the perimeter of the basin 'p'.

R1 = (𝐴𝐿𝑇𝑀𝑎𝑥-𝐴𝐿𝑇𝑀𝑖𝑛) 𝑃 = 𝑟 𝑝 (3)
Where ALTMax is the altitude of the culminating point and ALTMin is the altitude of the outlet. 3 -R2 is the dimensionless Schumm relief ratio [START_REF] Schumm | Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey[END_REF][START_REF] Melton | An analysis of the relations among elements of climate, surface properties, and geomorphology[END_REF][START_REF] Rawat | Hydro-morphometrical analyses of sub-himalyan region in relation to small hydro-electric power[END_REF] defined as the ratio between the total relief of a basin 'r' and the longest dimension of the basin parallel to the principal drainage line 'L':

R2 = 𝑟 𝐿 (4)
'L' is defined as the length of the equivalent rectangle:

L = 𝐾 𝐺 √𝐴 1,12 [1 + √1 -( 1,12 𝐾 𝐺 ) 2 ] (5)
The basin relief ratio (R1 and R2) indicates the general slope of a drainage basin and the intensity of erosion process occurring on the slope of the basin [START_REF] Schumm | Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey[END_REF][START_REF] Asfaw | Quantitative analysis of morphometry on Ribb and Gumara watersheds: Implications for soil and water conservation[END_REF]. 4 -The last topographic parameter is the Ruggedness number 'H' calculated as the product of drainage density 'Dd' and relief of basin 'r' both in the same units. It refers to the level of smoothness and roughness of the basin surface and its vulnerability for soil erosion. High ruggedness values indicate areas with uneven topography and sensitive to soil erosion [START_REF] Strahler | Statistical Analysis in Geomorphic Research[END_REF][START_REF] Asfaw | Quantitative analysis of morphometry on Ribb and Gumara watersheds: Implications for soil and water conservation[END_REF].

H = r Dd (6)
5 -'Mean-slope' (Fig. 2b) denotes the average slope of each basin, and was calculated from the DEM by the slope tool available on the GIS software. The procedure calculates for each cell the maximum rate of variation of the altitude values with respect to its neighbors [START_REF] Burrough | Principles of Geographical Information Systems[END_REF]. 6 -The Langbein concavity index [START_REF] Langbein | Profiles of rivers of uniform discharge[END_REF] was calculated as the ratio of the difference between ALTMax and H, the altitude at 50% distance between the source and the outlet, and the relief of a basin 'r'. The concavity index 'CI' indicates the channel steepness, the higher the value, the faster the channel gradient decreases. CI =

𝐴𝐿𝑇𝑀𝑎𝑥-𝐻 𝑟

(7) 2.2.4. Drainage density parameter: The drainage density is a key parameter that relates the operating geomorphological processes to the topography of a drainage basin [START_REF] Sangireddy | Controls of climate, topography, vegetation, and lithology on drainage density extracted from high resolution topography data[END_REF]. It is expressed as the ratio of total channel length to total basin area (Km/Km²). After extracting the drainage density, the average drainage density 'MEAN-Ddrain' (Fig. 2c), which represents the average drainage density in each basin, was extracted using the 'Area statistics' tool available on GIS software that calculates the values of the drainage density raster within each basin.

Geological parameters

The basin geological characteristics are represented by four parameters. First, the 'Structural control' parameter, which reflects the type of geological structuring that controls the drainage network and the basin geomorphology. Four types of structural controls were identified in the study area, by visual interpretation and image processing of satellite images, geological maps and other USGS data [START_REF] Van Der Meer | Potential of ESA's Sentinel-2 for geological applications[END_REF]Adiri et al., 2017;[START_REF] Fal | The contribution of Sentinel-2 satellite images for geological mapping in the south of Tafilalet Basin (Eastern Anti-Atlas, Morocco)[END_REF][START_REF] El Janati | Application of remotely sensed ASTER data in detecting alteration hosting Cu, Ag and Au bearing mineralized zones in Taghdout area, Central Anti-Atlas of Morocco[END_REF], but also from previous work [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF]. These four types of structural control are Tabular = 0, Folding = 1, Faulting = 2, and Faulting + folding = 3. The second parameter of the geological characteristics reflects the lithologic homogeneity or heterogeneity of the basin area. If a basin has a surface 80% covered by a lithological formation of the same lithological formation, it is described as homogeneous and it is assigned the code 0, otherwise, it is qualified as heterogeneous and it is assigned the code 1. The rock type parameter 'Rock-type' was based on a simplified basin lithology. For each basin, rock types were lumped into four major group, to which codes have been assigned: granitic-metasedimentary=3, volcanic=2, sedimentary=1, and quaternary sediment=0, broadly reflecting a change from harder to softer rock types. The parameter "Lineament density" reflects the average density of geologic lineaments identified in each basin, expressed as the ratio of total lineaments length to total catchment area (Km/Km²) (Fig. 2d). The structural alignments identification was performed automatically by a computer algorithm based on the Hough transform method [START_REF] Karnieli | Automatic Extraction and Evaluation of Geological Linear Features from Digital Remote Sensing Data Using a Hough Transform[END_REF][START_REF] Poncelet | Transformée de Hough et détection de linéaments sur images satellitaires et modèles numériques de terrain[END_REF][START_REF] Aydogan | Extraction of lineaments from gravity anomaly maps using the gradient calculation: Application to Central Anatolia[END_REF][START_REF] Renson | Automatisation de la détection des cratères lunaires sur des images et MNT planétaires[END_REF][START_REF] Bouramtane | Automatic Detection and Evaluation of Geological linear Features from Remote Sensing Data Using the Hough Transform Algorithm in Eastern Anti-Atlas (Morocco)[END_REF]. The Hough transform algorithm was applied to the multi spectral image Landsat 8 OLI satellite (30 m resolution), the image of the Sentinel-2A satellite band 11 (10 m resolution) and the radar image Sentinel-1A band Intensity-VH. The 3 sets of geological lineaments obtained on the 3 images were merged while eliminating the duplicate or triplicate to avoid repetition of the information. The relationships between some parameters taken 2 to 2 are presented in Fig. 5 that shows, for example, that the relief ratio 'R1' increases with the ruggedness number 'H' on all basin types, and that the scatter plots overlap widely. The parallel-type networks are slightly individualized (Fig. 5a andc), but their superposition with angular and lattice types networks does not allow them to be separated into an isolated group, as shown in the means and standard deviations cross plot (Fig. 5c). The same is true, for example, for the 'Lineament Density' and 'Mean Slope' parameters (Fig. 5b andd). This approach is therefore unsuccessful for the discrimination of types of drainage networks, hence the use of multivariate methods. The results of the PCA are presented in Fig. 6 and Table 2. The variance distribution according to the Principal Components (PCs) showed that the information was distributed in many PCs, indicating complex and nonredundant information. The first five PCs showed eigenvalues higher than unity indicating that they convey more information than a single variable. The first four PCs accounted for 72% of the total variance of the dataset. The first PC (PC1), which explains 35.5% of the total variance, essentially represented topographic parameters (H, R1, R2 and ALT-MAX), slope parameters (Mean-slope) as well as the parameter lineament density. Although less marked, the parameters related to the basin shape KH, the Homogeneity, the area, structural control and the type of drainage patterns were positively correlated to this PC1. The second PC, PC2, accounted for 14.4%, involving a group of morphometric parameters (Surface, Perimeter, KG) and to a lesser extent a second group of parameters related to type of drainage network, and the geological parameters structural control and Homogeneity. The type of drainage pattern was clearly correlated with the third PC (PC3, 12.7% of the variance) as well as the geological parameters (structural control and homogeneity) reflecting a close relationship between the type of drainage network and the geological parameters. This relationship seemed to be stronger for small and medium size basins as suggested by the negative correlation with the morphometric parameters (Surface, Perimetre and Concavity). Finally, the fourth PC (PC4, 9.3%) was correlated with the drainage density, secondarily with the morphometric parameter KG, the Area and the Ruggedness number H. Thus, the variance distribution showed a complex pattern, but the type of drainage pattern appeared on the first three PCs with strong correlation with topographic and geological parameters. For the latter, the contribution of the structural control and homogeneity played a major role. 4 shows the correlation values between the different parameters and the three discriminating functions. The first discriminant function (F1, 69.18%) relied mainly on the geological characteristics, and presented high correlation values for the parameters "structural control" and "homogeneity" (0.844, 0.933 respectively). The second discriminant function (F2, 24.20%) have a strong negative correlation with topographic parameters H, R1, R2 and ALT-MAX (-0.737, -0.588, -0.481 and -0.502, respectively) and with a less extend the drainage density parameter (Mean-Ddrain, -0.426). The third discriminant function (F3, 6.62%) relied mainly on the geological parameter "Area" with a correlation of 0.657, secondarily on the parameter that describe the basin mean slope (Mean-Slope, 0.576), and on the geological parameter "lineaments density" with a correlation of 0.515. Full equations of the discriminant functions are given in supplementary material S3. F1 discriminated very well the dendritic from the parallel pattern and from the group trellis + angular pattern. F2 separated the parallel pattern from the other patterns. Although F3 carried the least variance, it was the only function to separate trellis from angular patterns. The F1-F2 score plot (Fig. 7) displays each group with an accuracy rate of 93.6%. The discriminant function F1 separates the basins controlled by the geological factors with the structural control and homogeneity parameters, i.e. the angular and trellis patterns, from the basins not controlled by this factor, i.e. the dendritic and the parallels. The centroids of angular and trellis patterns have positive coordinates along F1 axis, whereas those of dendritic and parallel pattern are negative and close to zero, respectively (Fig. 7). F2 separates the basins that are influenced by a high topography and high relief ratio, as is the case with the parallel pattern, from the basins that are weakly influenced by this factor, for example the dendritic pattern, or those not controlled by this factor, such as trellis and angular patterns [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF].

The F1-F3 plot (Fig. 7)

DISCUSSION

Trellis and angular drainage patterns are widespread in the eastern Anti-Atlas, suggesting that hydrological processes in this region are subject to structural control [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF][START_REF] Stokes | Active and passive tectonic controls for transverse drainage and river gorge development in a collisional mountain belt (Dades Gorges, High Atlas Mountains, Morocco)[END_REF][START_REF] Pereira-Claren | Planform geometry and relief characterization of drainage networks in high-relief environments: An analysis of Chilean Andean basins[END_REF]. The first 3 principal components highlight the parameters that describe the topography, slope and geological parameters (Structural control, Area, Homogeneity and Lineaments density), all positively correlated with the type of drainage pattern. These results are in agreement with the work of [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF], [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF] and [START_REF] Twidale | River patterns and their meaning[END_REF], as the factors that explain the most variance between the different types of drainage pattern are those that reflect the topography and geological structuring [START_REF] Stokes | Active and passive tectonic controls for transverse drainage and river gorge development in a collisional mountain belt (Dades Gorges, High Atlas Mountains, Morocco)[END_REF][START_REF] Pereira-Claren | Planform geometry and relief characterization of drainage networks in high-relief environments: An analysis of Chilean Andean basins[END_REF].

Principal components furnish macro parameters, i.e. synthetic data that convey strong and significant information and, therefore, are particularly suitable for digital mapping and spatial analysis [START_REF] Rezende-Filho | Impacts of lithological and anthropogenic factors affecting water chemistry in the upper Paraguay River Basin[END_REF]. The distribution of PCs values is showed in Fig. 10. The Highest values of PC1 denote the influence of a high relief with steep slope, heterogeneous rock lithology and high intensity of fracturing and geological deformation. Higher PC1 values are observed in the SOA sector, on the one hand, at the heart of the Neoproterozoic inliers where angular drainage patterns occur, and on the other hand on the northern and southern flanks of the Saghro inlier, where strong slope gradient are associated with parallel patterns. Moderate values along PC1 are found in the MB, OOA and BS sectors, where almost parallel or slightly arcuate outcrops exert structural control on the drainage pattern and reduce the influence of slope parameter. This led to the formation of trellis patterns associated with the Bani structures in the BS and MB sectors, and with anticlinal and synclinal flanks in folded sedimentary sequences in the OOA sector. Low and negative values along PC1 reflect low relief, a relatively gentle slope and a lithology characterized by Cretaceous and Quaternary sedimentary deposits, and a lack of prominent geologic and tectonic structures. They are observed in the KKD and STB sectors, reflecting the weak influence of the slope factor and the presence of dendritic models. The vast alluvial zones of these areas strongly influence this phenomenon. The observation of some parallel network in the western KKD sector indicates a slight change in slope gradient from east to west. [START_REF] Jung | Under what conditions do parallel river networks occur?[END_REF] observed a similar phenomenon and calculated that the drainage networks present a transition from dendritic to parallel from an increase of 1% to 3% of initial regional slope. PC2 is also associated with the type of drainage network. It is positively correlated with morphometric parameters (Surface, perimeter and KG) and geological parameters (Structural Control and Homogeneity) and negatively correlated with drainage density. These characteristics depend mainly on the size and shape of the basin, the lithology and the extent of the drainage network [START_REF] Horton | Remarks on hydrologic terminology[END_REF][START_REF] Melton | An analysis of the relations among elements of climate, surface properties, and geomorphology[END_REF]. High PC2 values (> 4) are observed at northeastern and western part of the KDD sector, in the northwestern and southwestern parts of the Erfoud region, east of Ikniouen, and on the Bani structures in the east -northeast of Zagora (Fig. 10), where most of the basins are large and elongated. Negative values on PC2, which reflect a dependence on drainage density, were identified in the south and north flanks of the SOA sector and in the southwestern part of the BS sector. In these areas, steep slopes on clay, shale and sandstone rocks result in higher drainage density [START_REF] Melton | An analysis of the relations among elements of climate, surface properties, and geomorphology[END_REF]. In the STB sector south of Erfoud, the negative values correspond to the alluvial plain, where the absence of geological structuring and low slope lead to high drainage densities.

PC3 highlights the relationship between the geological parameters and the type of drainage pattern. These parameters collectively reflect the influence of the geological structuring on the control of drainage network for the basins with the angular and trellis pattern. In other words, the basins with a high coordinate in PC3 are likely to reflect the presence of one of these two types of drainage patterns. The high positive PC3 values are observed in the MB, OOA, and in the east of BS, SOA and STB sectors, where both trellis and angular patterns prevail on small to medium size basins with relatively heterogeneous lithology, over folded and fractured terrain (Fig. 10). Negative values in PC3 reveal a dependence on topographical characteristics, in particular altitude variation, and basin size. Such negative values occur in KKD and east STB sectors in wide basins with gentle slope and mainly dendritic networks, on Cretaceous deposits in the first and alluvial plains in STB sector, attesting for the weak influence of structural control on the drainage network. Generally large basins with branching networks like apple trees indicate a lack of control structure or slope over the drainage pattern, which is observed in the dendritic pattern [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF]. The low PC3 values (between 2 and 0) in the western part of SOA sector is related to the high variation in altitude from the center of the Saghro inlier towards its northern flank. In this region, hybrid drainage patterns occur, from angular networks in wide basins upstream, following major faults, to parallel networks with narrow basins downstream in the flank of the inlier. [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF] observed the same behavior in northwestern Finland, where some rivers run parallel to fault zones. These faults are structural and often accompanied by conjugated fractures, leading an angular configuration of the drainage network.

Finally, high PC4 values occur in areas that presents a combination of high drainage density and the presence of sedimentary deposits, which is particularly clear in the southern flank of the SOA sector and the northern part of the Zagora region. Together, the results highlighted by PC2 and PC4 are consistent with observations made by [START_REF] Day | Lithologic controls of drainage density: A study of six small rural catchments in New England[END_REF], [START_REF] Mueller | Sediment supply and channel morphology in mountain river systems: 1. Relative importance of lithology, topography, and climate[END_REF] and [START_REF] Sangireddy | Controls of climate, topography, vegetation, and lithology on drainage density extracted from high resolution topography data[END_REF], as sedimentary rocks show relatively high values of drainage density compared with granitic, volcanic, and meta-sedimentary rocks.

The three discriminant functions Identified from the Linear Discriminant Analysis show that the basins can be distinguished with a high rate (90%) from their type of drainage pattern, using geologic, topographic and morphometric parameters. The confusing patterns mainly involve trellis and angular patterns and are false positive classifications. Indeed 12 trellis patterns were discriminated as angular and 8 angular patterns were discriminated as trellis. In the same way F1, strongly correlated with the parameters that describe the geological structuring, groups the trellis and angular patterns in a same set and distinguishes them from the dendritic and parallel patterns, characterized by the absence of geological structuring. Both analyses, PCA and LDA, provide consistent results, especially PC3 and F1, given the close relationship between trellis and angular patterns, and the similar geological and geomorphological conditions for their formation. These conditions are folded sedimentary branching, converging, diverging or roughly parallel faults and joints and/or faults at right angles, which are generated by similar geological structuring factors [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF][START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF][START_REF] Twidale | River patterns and their meaning[END_REF]. Finally angular and trellis patterns were separated by F3, because of the high weight of the Area parameter, representing the geological sector of the basin's location. The distinct geology and geomorphology of each sector play a major role in determining the type of drainage systems, particularly those heavily influenced by geological structuring. Angular pattern basins are mainly located in the Saghro-Ougnat sector, i.e. areas composed of Proterozoic and Paleozoic terrains with high fracturing densities. This area was strongly impacted by several phases, mainly during the Neoproterozoic Panafrican and the late Paleozoic Hercynian orogenies remobilizations [START_REF] Hefferan | A reconsideration of Pan African orogenic cycle in the Anti-Atlas Mountains, Morocco[END_REF]. Many of the inherited faults were reactivated during the Alpine collision and High Atlas inversion, this last event did not cause any direct deformation in the Eastern Anti-Atlas, except a general uplift combined with the effect of a recent thermal anomaly, these inherited faults still constituting weak features [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Frizon De Lamotte | The Atlas System BT -Continental Evolution: The Geology of Morocco: Structure, Stratigraphy, and Tectonics of the Africa-Atlantic-Mediterranean Triple Junction[END_REF]Soulaimani et al., 2014;[START_REF] Gouiza | Post-Variscan evolution of the Anti-Atlas belt of Morocco constrained from low-temperature geochronology[END_REF]. Among these faults in SOA sector, Bou laghzazil-Tinifift, Tizi n'boujou and Arekouz fault (Fig. 11a andb) are the main ones controlling a drainage network that follow the orientation of faults corridors, drawing angular shapes. Our finding are consistent with the conclusions of [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF], as angular drainage networks occur in area with prominent structural control, and the pattern is directly conditioned by the right-angled jointing or faulting of rocks [START_REF] Stokes | Active and passive tectonic controls for transverse drainage and river gorge development in a collisional mountain belt (Dades Gorges, High Atlas Mountains, Morocco)[END_REF][START_REF] Burr | Morphology of fluvial networks on Titan: Evidence for structural control[END_REF][START_REF] Pereira-Claren | Planform geometry and relief characterization of drainage networks in high-relief environments: An analysis of Chilean Andean basins[END_REF].

On the other hand, trellis pattern networks are mainly found in the Bani (BS), Maider basin (MB) and Ougnat-Ouzina axis (OOA) sector (Fig. 11d-c), i.e., in areas with a high deformation of geological structures and the presence of huge folded outcrops. Trellis patterns occur around and on the flanks of the prominent folds in the OOA sector, some with spectacular shape consisting in large crescent-or boomerang-shaped Cambrian structures, easily observed on satellite images. They are the Jbel Tijekht, Jbel Tadout and Shayeb Arass anticlines, as well as Ottara and Amesseoui synclines. [START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF] reported that these large Cambrian crescent-shaped or boomerang structures reflect inversion tectonics, describing locally a double fold mechanism controlled by faults that lead to basement exposure. Similarities and classification confusion between angular and trellis patterns has been observed by other authors. Based on the deviations of drainage patterns from self-similarity, [START_REF] Mejía | Identification and characterization of dendritic, parallel, pinnate, rectangular, and trellis networks based on deviations from planform self-similarity[END_REF] observed that the shape characteristics of dendritic networks are self-similar at different scales of observation. The shape characteristics of parallel networks are self-affine, presenting an invariant image by anisotropic transformation and both angular and trellis networks are approximately self-similar. More recently, Jung et al. ( 2019) used a classification method based on quantitative parameters, namely the tributary junction angles, discriminate 5 drainage network patterns classes, dendritic, parallel, pinnate, trellis and rectangular. Their proposed method clearly discriminated dendritic, parallel, and pinnate networks, but differences between angular and trellis network types were hardly distinguished.

In the Banis sector, the trellis patterns occur in a geomorphological context that consists of alternating cuestas and deep depressions, typical of the EAA. These cuestas consist of sandstone ridges of the first and second Bani of Ordovician-Silurian age. The drainage networks flow in a straight line along to the greatest slope on the backs of the Cuestas towards the Feijas depressions. Then, the drainage networks adopt their layout and their orientation, which gives a characteristic trellis shape, controlled by the edges of the outcrops of formations, whose resistance varies in the parallel belts [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF] (Fig. 11d). The result is consistent with the idea that the topography is the determining factor of this type of drainage pattern, which is generally associated with steeply to moderate sloping regions [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF][START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF][START_REF] Twidale | River patterns and their meaning[END_REF][START_REF] Radaideh | Detection and analysis of morphotectonic features utilizing satellite remote sensing and GIS: An example in SW Jordan[END_REF]. In the study area, parallel pattern basins with relatively high relief values (Fig. 2a) are located in the northwest, on the southern edge of the Saghro inlier, and to the southwest on the escarpment of the Feijas, and northeast of the city of Zagora (Fig. 12a andb). They are also found west of the KKD sector, with higher maximum altitude and slope gradient (Fig. 2a andb). In addition, the parameters that reflect the drainage density (MEAN-Ddrain) have high coordinate values in F2. This suggests the existence of a process that, in addition to the control by the slope, incorporates the drainage density. This observation is consistent with [START_REF] Johnson | Streams and Their Significance[END_REF], who suspected a relationship between slope and drainage density: the steeper the slope, the more parallel patterns are developed, and the more the drainage density increases [START_REF] Twidale | River patterns and their meaning[END_REF], as it appears in Figure 11 (a and b). The discriminant function F2 makes it possible to distinguish basins as a function, on the one hand, of the influence of the slope on their drainage networks and, on the other hand, of the drainage density. Thus, F2 makes it possible to separate the parallel pattern basins controlled by high relief and a steep slope gradient, the dendritic pattern basins characterized by a low slope gradient, and finally the trellis and angular pattern basins, mainly controlled by the geological structuring.

A few confusions in the parallel patterns were made using discriminant functions, with two parallel patterns identified as trellis, and one angular as parallel. In the SOA sector, angular networks are observed within the inliers in a strongly fractured area, and they connect to parallel networks developed on the strongly inclined Palaeozoic cover in the southern edge of the inliers. Such organization of the drainage networks may reflect a structural control by the formation of the Inliers. The relationship between angular and parallel in region where structural control occur has already been described by [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF] that observed rivers running along parallel fault in north-western Finland. These faults, structural and most commonly accompanied by conjugate fractures, lead to angular network pattern. In the Bani sector, parallel and trellis networks developed on the same folded Palaeozoic formations, however they differ in their degree of parallelism resulting from a progressive steepening of slope, in relation with the degree of resistivity to erosion (Fig. 11d and 11b) [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. In the OOA sector, trellis drainage networks are restricted to a parallel relationship between tributaries (Fig. 11c). In the simple parallel model, parallelism may occur between tributaries or unrelated streams (Fig. 12a andb) [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF]. As revealed by PC3, in the Bani sector, because of alternation between ridges of hard sandstone (first and second Bani) and clay levels (Feija and Ktaoua shales) trellis networks have developed, although restricted to a parallel relationship between secondary tributaries, which are usually elongated and connect downstream at right angles [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF] (Fig. 11c). In addition, parallel models have developed north of Zagora on the escarpment of the first Bani, flowing southwards over the gently inclined clay deposits. Because of uniform resistivity to erosion, the drainage network show a parallelism between all tributaries as well as between unrelated channels (Fig. 12b).

The separation of the dendritic patterns from others patterns in the two-dimensional space F1-F2 reflects a specificity of this drainage pattern. The positioning of the centroid with a positive but close to zero coordinate on F2, and clearly negative on F1 reflects the control by a low slope gradient and the absence of control by geological structuring. Of course, dendritic-pattern basins were not expected to have high F1 values, because these drainage patterns are formed in areas without significant structural control, a lack of geomorphologic structures, a homogeneous lithology, and a slope gradient from weak to zero [START_REF] Zernitz | Drainage Patterns and Their Significance[END_REF][START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. These characteristics generate dendritic patterns characterized by irregular branching in all directions with the tributaries joining the main stream at any angles. These findings are observed in the EAA, where the dendritic patterns are exclusively in tectonically stable zones, with low relief and homogeneous lithology, as the case of the Kem-kem domain, consisting exclusively of Cretaceous deposits, and within the sedimentary basin of South-Tafilalat with flat lying Quaternary sediments (silt, sand, alluvial fans and talus) (Fig. 12c andd). Due to lack of marked structural control, and homogeneous lithology in these areas, the drainage network consists of channels oriented in a wide range of directions, leading to high drainage density. The results of the Agglomerative Hierarchical Clustering analysis are in agreement with the results of the LDA with respect to, on the one hand, the confusion between the parallel, trellis and angular patterns, due to similar factors responsible for their formation, and on the other hand, the individualization of a dendritic pattern. The position of the parallel type in the dendrogram between the angular+trellis group and the dendritic patterns could not be only due to its formation conditions, but could also stem from transition phenomenon from one drainage pattern to another over time or space, as mentioned by [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. This author asserts that all transitions are possible from parallel to trellis or dendritic patterns and vice versa.

Here, the main transition in drainage patterns in space occur in the north-south or northeast-southwest direction, from angular pattern in the Saghro and Ougnat inliers (Fig. 11a andb) to trellis patterns in the Bani, MB and OOA sectors (Fig. 11c andd), and then to dendritic patterns in the KKD sector (Fig. 11d). Such transition corresponds to the passage from the Anti-Atlasic domain represented by the SOA sector to the Ougartian domain related to the Ougarta chain. The Ougartian domain interferes with the anti-atlas in the OOA sector, continues under Cretaceous and Neogene formations in the KKD sector, and ends abruptly in the north, crosscut by the Oumjerane-Taouz Fault. It outcrops further in the southeast [START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF]. The transition in time from parallel to trellis pattern can occur during a progressive structural control due to the development of erosion [START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. In the western part of the KKD sector, the drainage network is shifting from dendritic to parallel. This may be due to the slow uprising and increase in the slope gradient that have occurred in the EAA since the Lower Mesozoic to the present day [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF][START_REF] Jung | Under what conditions do parallel river networks occur?[END_REF][START_REF] Soulaimani | Late Ediacaran-Cambrian structures and their reactivation during the Variscan and Alpine cycles in the Anti-Atlas (Morocco)[END_REF][START_REF] Gouiza | Post-Variscan evolution of the Anti-Atlas belt of Morocco constrained from low-temperature geochronology[END_REF]. In this structural framework, changes in the orientation of the drainage system over geological time are indicators of transition from one event to another and reflect changes in the factors that control the formation of drainage patterns, as well as transition patterns [START_REF] Malik | Active tectonic influence on the evolution of drainage and landscape: Geomorphic signatures from frontal and hinterland areas along the Northwestern Himalaya, India[END_REF][START_REF] Brookfield | Evolution of the great river systems of southern Asia during the Cenozoic India-Asia collision: Rivers draining north from the Pamir syntaxis[END_REF][START_REF] Stokes | Active and passive tectonic controls for transverse drainage and river gorge development in a collisional mountain belt (Dades Gorges, High Atlas Mountains, Morocco)[END_REF][START_REF] Radaideh | Detection and analysis of morphotectonic features utilizing satellite remote sensing and GIS: An example in SW Jordan[END_REF]. The principal orientations of the drainage networks (Fig. 9) developed on Precambrian and Paleozoic (NE-SW, NW-SE and E-W) reflect the principal orientations of the geological structures over the entire EAA and mentioned in the literature [START_REF] Baidder | Fold interference pattern in thick-skinned tectonics; a case study from the external Variscan belt of Eastern Anti-Atlas, Morocco[END_REF][START_REF] Bouramtane | Automatic Detection and Evaluation of Geological linear Features from Remote Sensing Data Using the Hough Transform Algorithm in Eastern Anti-Atlas (Morocco)[END_REF]. They are inherited faults that operated during the collapse of the Pan-African belt and during Cambrian rifting. They have been reactivated as normal faults during the Paleozoic sedimentation and then reversed during the Hercynian event. At the end of the Hercynian event, most of the structure of the Anti-Atlas was established, apart from the reactivation of the inherited faults [START_REF] Soulaimani | Late Ediacaran-Cambrian structures and their reactivation during the Variscan and Alpine cycles in the Anti-Atlas (Morocco)[END_REF]. During the Lower Mesozoic and Lower Neogene, the Eastern Anti-Atlas did not record any significant deformation [START_REF]Géologie structurale de l'Anti-Atlas oriental[END_REF]. However, this period was characterized by a succession of subsidence and uplift events that gave the current topography of the EAA, in particular with higher altitude in the Saghro and Ougnat massifs [START_REF] Gouiza | Post-Variscan evolution of the Anti-Atlas belt of Morocco constrained from low-temperature geochronology[END_REF] and a gradual decrease towards the south and southeast. The directions observed in the rose diagrams for networks developed on the most recent formations rather reflect the influence of the slope gradient. These observations suggest that the drainage system shifts from the influence of structural control in Precambrian and Palaeozoic formations to the influence of topography in the post-Hercynian events from the Lower Mesozoic to the Quaternary. Such observations are consistent with the spatial and temporal transition of the drainage patterns.

CONCLUSION

In this work, we have studied the distribution and arrangement of the drainage network in the Eastern Anti-Atlas, an area where the aridity facilitates satellite images of surface formations and drainage networks. We assumed that the patterns of drainage network were not randomly arranged and we looked for factors that could control them. A dataset of 230 basins and 16 parameters related to the topography, morphometry, geology was analyzed by conventional multivariate statistical methods such as principal component analysis (PCA), linear discriminant analysis (LDA) and agglomerative hierarchical clustering (AHC). The number of basins studied is large enough to achieve a reliable statistical approach. The PCA shows that the 16 parameters can be summed up in a rather limited number of independent sources of variability, since the first 6 PCs carry about 86.06% of the information. The number of independent sources of variability ( 6) is small compared to the number of objects studied (230). Four types of drainage network patterns have been identified, namely trellis, angular, parallel and dendritic patterns. The information contained in this dataset is complex and not redundant. The variation between the different drainage patterns is first explained by parameters related to the topography, the slope and the geology. The LDA made it possible to distinguish between the four types of drainage patterns with a success rate of 90%, using 3 discriminant functions. The dendritic pattern is associated with tabular and homogeneous sedimentary geological formations, little fractured with a small slope gradient. The trellis and angular patterns are associated with geological structures that exhibit intense deformation and fracturing in the study area. They are discriminated based on the distinct geological and geomorphological characteristics of each geological sector. Finally, the parallel pattern network is associated with steeper slopes. Moreover, the similar results from both, the linear discriminant analysis and the unsupervised hierarchical clustering, show confusion and similarity between the parallel, trellis and angular patterns, due to, on the one hand, the similar geological setting responsible for their formation and on the other hand, the possible transition (in time and space) between these three patterns in the study area. It suggests a possible spatial and temporal transition between the four types of drainage patterns, being the result of a possible shifting of the drainage network, from the influence of structural control in Precambrian and Paleozoic formations, to the influence of topography from post-Hercynian events in the Lower Mesozoic to the Quaternary. Principal Component Analysis, Linear Discriminant Analysis and agglomerative hierarchical clustering are reference statistical methods, used widely in science. However, the approaches used here suffer from a limitation, in the sense that the data used to search for discriminant functions were calculated on the entire dataset, i.e. on all 230 basins. The discrimination success rate, that is, the validation of the performance of the linear discriminant analysis, is measured on the same set of data using the confusion matrix. A methodological study based on the separation of the data into, on the one hand, a set for the determination of learning methods of the type of drainage patterns from the different parameters, and on the other hand, an independent validation set, would be more in line with a machine learning-type approach. A methodological study to compare the results of these two approaches is currently being carried out and should be presented shortly.

  Statistical approach of factors controlling drainage network patterns in arid areas. Application to the Eastern Anti Atlas (Morocco) Tarik Bouramtane, Suzanne Yameogo, Meryem Touzani, Abdessamad Tiouiouine, M'Hamed El Janati, Jamila Ouardi, Ilias Kacimi, Vincent Vallès, Laurent Barbiero To cite this version: Tarik Bouramtane, Suzanne Yameogo, Meryem Touzani, Abdessamad Tiouiouine, M'Hamed El Janati, et al.. Statistical approach of factors controlling drainage network patterns in arid areas. Application to the Eastern Anti Atlas (Morocco). Journal of African Earth Sciences, 2020, 162, pp.103707. 10.1016/j.jafrearsci.2019.103707. hal-02389502

  Fig. 2: Map of studied basin in the eastern Anti-Atlas, Green dots indicate locations of main cities, basins are outlined in bold: (a) Elevation (b) Mean basin Slope (c) drainage network mean density and (d) geological lineaments mean density.

Figure 3 :

 3 Figure 3: Distribution of the different drainage network patterns in the Eastern Anti-Atlas. (See figure 4 for black boxes a, b, c and d) 3. RESULTS Drainage network distribution: A total of 230 basins with area ranging from 9.5 Km² to 482 Km² were extracted. The drainage network patterns were identified as follows (Fig. 3). One hundred and six basins showed trellis drainage patterns developing mainly on the Paleozoic folded cover. Fifty basins had angular patterns and developed on the ancient Precambrian formations in both Saghro and Ougnat inliers, and in the extreme northwest on the Erfoud anticlinorium. Forty five basins were identified as dendritic patterns in the Ziz alluvial plain to the northwest, and in the southwestern extent of the Kem-Kem domain. Finally, twenty nine basins had parallel

Fig. 5

 5 Fig. 5: a) and b) : Distribution plot for the parameters Rudgeness number 'H' and Relief ratio 'R1', and Mean-Slope and Lineament Density, respectively, for the 230 extracted basins and depending on the type of drainage network. c) and d) Mean values and standard deviations for the same parameters.

  Fig. 7: Distribution of the observations and centroids of each type of drainage network pattern in F1-F2 and F1-F3 score plots (0 = Dendritic, 1 = Parallel, 2 = Trellis, 3 = Angular). Agglomerative Hierarchical Clustering: PCs are a linear combination of the 16 variables, i.e. macro-parameters that can be used for classification. The first seven PCs accounted for about 90% of the variance of the dataset. The results of the agglomerative hierarchical clustering based on the coordinates of the centroids throughout the axes of the first seven PCs are presented in Figure8. This made it possible to specify the degree of similarity or dissimilarity between the different types of drainage network patterns. Clearly, the dendritic pattern was individualized from the other patterns, whereas the trellis and angular patterns showed higher similarity, forming an isolated group. The parallel pattern was halfway between the group angular + trellis and the dendritic patterns.

Fig. 9 :

 9 Fig. 9: Change in the azimuth direction of the drainage network in the Eastern Anti-Atlas from Precambrian to Quaternary formations.

Fig. 11 :

 11 Fig. 11: Examples of relationship between major structural features (faults, synclines, anticlines) and drainage networks in (a and b) the Saghro-Ougna axis sector, (c) the Ougnat-Ouzina axis sector, and (d) the Bani sector. AOF = Assif n'Oussif Fault; BTF = Boulghzazil-Tinifift Fault; AF = Akrouz Fault; TRF = Tizi n'Ressas Fault; ASA = Anticline of Shayeb Arras; SA = Syncline of Amessoui ; AT = Anticline of Tijekht; B1 and B2 = Sandstone first and second Banis. (See the locations of the maps in Supplementary Material S4)

Fig. 12 :

 12 Fig. 12: Examples of relationship between parallel and dendritic drainage patterns, and slope gradient, geological structures and drainage density. (a) Parallel drainage patterns on the western edge of Saghro inlier in Saghro-Ougnat axis; (b) Parallel drainage network north of Zagora city in the Feija escarpment (BS Sector). (c) Dendritic drainage patterns in the Kem-Kem Domain, (d) Dendritic drainage patterns in the South-Tafilalt basin. (See the locations of the maps in Supplementary Material S4)The results of the Agglomerative Hierarchical Clustering analysis are in agreement with the results of the LDA with respect to, on the one hand, the confusion between the parallel, trellis and angular patterns, due to similar factors responsible for their formation, and on the other hand, the individualization of a dendritic pattern. The position of the parallel type in the dendrogram between the angular+trellis group and the dendritic patterns could not be only due to its formation conditions, but could also stem from transition phenomenon from one drainage pattern to another over time or space, as mentioned by[START_REF] Howard | Drainage Analysis in Geologic Interpretation: A Summation[END_REF]. This author asserts that all transitions are possible from parallel to trellis or dendritic patterns and vice versa. Here, the main transition in drainage patterns in space occur in the north-south or northeast-southwest

Table 1 :

 1 Summary details of the 6 sectors

	Sector	Area	Mean	Mean	Dominant lithology Geomorphological
		(Km²)	Elevation (m)	Slope (°)		Characteristics
	BS	7810	882	7,9	Sandstone and quartzite, micaceous	First and second Bani ridge: monoclinal structures,
					sandstone clay, green	Tazzarine basin
					shale
	OOA	3527	780	7,7	Green sandstone and quartzite, green	shear zone that consists of a succession of anticlines and
					sandstone and volcanic	synclines
					rock lenses, quartzite
					sandstone, green shale
					siliceous sandstone,
					marl limestone,
					sandstone and pelite
	MB	3120	780	6,2	Sandstone and quartzite, schist, pelitic	basin structures and depressions, such as the
					schist	Fezou syncline, lowland area
						at the foothills of the high
						reliefs forming part of a
						sedimentary basin
	KKD	5048	743	6,5	Limestone, Marl, Siliceous limestone	Tabular sector of Cretaceous deposit
	STB	2717	795	5,4	pelite, sandstone and quartzite, limestone,	lowland area, plains, mostly synclinal in character,
					alluvial conglomerate,	anticlinal structures in the
					coarse sand, limestone	northern and southern part
					cement conglomerate	of the area
	SOA	9255	1412	12,2	Amphibole granite, Rhyolite, Basalt, shale,
					sandstone,
					conglomerate, pink
					sandstone, green shale.

Table 2 :

 2 Principal Components, respective eigenvalues and contribution to the total variance.

		PC1	PC2	PC3	PC4	PC5	PC6
	Eigenvalues	6.04	2.45	2.16	1.59	1.35	0.94
	Contribution (%)	35.53	14.42	12.69	9.34	7.96	5.50
	Cumulated contribution (%)	35.53	49.94	62.63	71.97	79.93	85.44

Table 4 :

 4 Correlations between variables and discriminant functions (higher values are in bold)

		F1	F2	F3
		(69.2%)	(24.2%)	(6.6%)
	H	0.288	-0.737	0.267
	R1	0.294	-0.588	0.293
	R2	0.299	-0.481	0.250
	Alt max	0.343	-0.502	0.379
	Concavity	-0.326	0.354	0.357
	Area	0.222	-0.281	0.657
	Structural control	0.844	-0.007	0.147
	Homogeneity	0.933	0.027	-0.033
	Rock type	0.284	0.181	0.056
	Surface	0.112	0.053	0.010
	Perimetre	-0.005	0.089	0.114
	KG	-0.231	0.120	0.154
	KH	0.196	-0.097	-0.098
	MEAN-Ddrain	-0.053	-0.426	-0.257
	MEAN-Slope	0.302	-0.176	0.576
	Lineaments density	0.344	-0.045	0.515


I Essaadi
email: imane.essaadi@emines.um6p.ma

B Grabot

P Féniès

Design of hybrid multimodal logistic hub network with postponement strategy

Keywords: Logistic hub, postponement, distribution chain, service level, hierarchical structure, multi-period, multi-commodity

This paper aims at suggesting a method allowing to design a logistic hub network in the context of postponement strategy, postponement being performed in hubs having industrial facilities in addition to logistic ones. We propose a two-stage mathematical mixed integer linear programming model for: 1) logistic hub network design 2) postponement location on the designed hub network. The suggested model manages characteristics not yet taken into account simultaneously in the literature: hierarchical logistic structure, postponement strategy, multi-commodity, multi-packaging of goods (raw materials or components vs. final products), multi-period planning. The solutions are compared through services levels and logistic costs.

Introduction

The implementation of networks of logistic hubs usually allows to decrease transportation costs and delivery delays in comparison with direct source/destination transportation [START_REF] Alumur | Network hub location problems: The state of the art[END_REF]. Within logistic hubs, material flows coming from different origins are sorted, consolidated depending on their destination then transported using unimodal or multimodal transport. Two families of hubs can be distinguished: pure logistic hubs, providing standard logistic services (warehousing, inventory management, packaging, labeling, orders preparation or cross-docking, sorting and transport /distribution) and combined logistic / industrial hubs, offering high addedvalue services on logistics (such as co-packing) and/or industrial functionalities allowing the final customization of the product. The late customization of products (often called "postponement") was embraced for many years by some industries such as computers, printers, medical products and fertilizers [START_REF] Wang | Case study on the application of postponement strategy and managerial insights[END_REF]. Pushing this logic to its limits, multinational firms now attempt to customize their products within their distribution centers, like Hewlett-Packard producing DeskJet printers in its factory in Singapore and customizes them for the European and Asian markets within its European distribution center near Stuttgart, Germany. In that new context, this paper aims at defining the optimal design of a network of logistic hubs with postponement strategy integration in a context of mass customization. The distribution network is supposed to be composed of four levels: production plants, regional logistic hubs, sub-regional logistic hubs, and urban/rural distribution centers. The network processes different goods with possibly different packaging (raw materials/components and finished products).

State of the art

The hub network design problem, also known as the distribution network design problem, was intensively studied in the literature on global distribution networks [START_REF] Alumur | Network hub location problems: The state of the art[END_REF][START_REF] Campbell | Twenty-five years of hub location research[END_REF][START_REF] Alumur | Hierarchical multimodal hub location problem with time-definite deliveries[END_REF]]. An analysis of this literature is summarized in Figure 1.

Figure 1. Logistic hub network review analysis

In the literature, few authors consider a multi-level network, like [START_REF] Alumur | Hierarchical multimodal hub location problem with time-definite deliveries[END_REF] while no paper addresses the integration of industrial services. Origin and destination nodes might be either allocated to a unique hub or to multiple ones: [START_REF] Gelareh | Multi-period hub location problems in transportation[END_REF] considered multiple allocation of clients to located hubs while [START_REF] Sohrabi | On Comparing Dedicated and Hyperconnected Distribution Systems: An Optimization-Based Approach[END_REF] studied multiple allocation of both plants and customers to intermediate hubs. Service level can be addressed through the definition of a maximum distance between distribution center and market zone [START_REF] Melo | Dynamic multi-commodity capacitated facility location: a mathematical modeling framework for strategic supply chain planning[END_REF] or delivery delay such [START_REF] Alumur | Multimodal hub location and hub network design[END_REF][START_REF] Rodriguez | Hub location under capacity constraints[END_REF]. Many papers assume the demand to be deterministic, which is seldom true, the demand usually evolving through time. Multiple commodity, allowing to consider products using different transportation means, is only considered by few researchers (cf. Figure 1). To our best knowledge, no paper considers a physical transformation of goods while transiting a hub (from bulk material to packs or pallets for instance). [START_REF] Albareda-Sambola | Multiperiod Location-Routing with Decoupled Time Scales[END_REF] state that one of the main characteristics of supply chain network design models is their multiperiod nature. Many studies have analyzed the advantages and disadvantages of various postponement strategies [START_REF] Yang | Implications of postponement for the supply chain[END_REF][START_REF] Van Hoek | The rediscovery of postponement a literature review and directions for research[END_REF][START_REF] Hsuan Mikkola | Supply-chain integration: implications for mass customization, modularization and postponement strategies[END_REF]. However, quantitative models for postponement implementation decisions are scarce: [START_REF] Schwartz | Designing distribution networks taking into account aspects of postponement, in The Development of Collaborative Production and Service Systems in Emergent Economies[END_REF] consider decisions on where to implement assembly and packaging functions in a distribution network. 2) while [START_REF] Schulze | Location-allocation model for logistics networks with implementing commonality and postponement strategies[END_REF] addressed the problem of facility location-allocation (plants, warehouses) considering commonality and postponement strategies in the logistic network. This quick analysis of the literature shows that no study gathers yet all the characteristics we have chosen to address in order to answer to present real problems, summarized by the "Proposal" column in Fig. 1.

Problem formulation

The design of a hybrid logistic hub network involving postponement strategy implies determining simultaneously the location of the logistic hubs and of the postponement services, while defining the optimal routing of flows minimizing total logistics costs. Postponement units will have as inputs raw materials and components coming from international plants, based on their specialization and logistic costs. They will provide bagged / assembled products in response to the requirements of the market zones. These requirements may differ in terms of packaging preferences and required response time. We assume that postponement units hold sufficient component inventories for meeting the customer deterministic demands. Furthermore, we do not make a priori assumptions on the capacity of the hubs and postponement units, since we consider a new design and not the reuse of existing facilities. The capacities (processing and storage) will be determined a posteriori (through simulation) by considering the expected levels of service. Market zones are allocated to a unique hub based on logistic costs. In this problem, the interdependence of the decisions makes it difficult to instantiate all the decision variables simultaneously. Indeed, the location of postponement units will impact the management of the logistic flow, as they constitute decoupling points. On the other hand, their location depends on the location of the logistic hubs and of the allocated demand (volume, response time and product preferences). For addressing this problem, we have chosen as a first approach to decouple the initial problem in two sub problems: definition of the logistic network, then location of the postponements facilities, even if the network structure may in theory be set into question by the positioning of the postponement units. Each sub-problem will be modeled by a deterministic mixed integer linear programming models (cf. §3.1 and §3.2)

First Sub-problem: logistic hub location problem

Within this sub-problem we have to decide: 1) the location of hub h among potential locations H using the zh binary variable. 2) The allocation of the origin and destination nodes to the located hubs, using transport mode m represented by the yo,d m binary variable. This latter is defined only if a modal link between node "o" and "d" exists i.e. Link (o, d, m) =1. 3) The flows routing within the network x m,k,t o,d i.e. the amount of final product k originated from plant p and transported from node "o" to node "d" using vehicle mode m and under packaging n at period t. The generated solution must provide the best benefit considering the initial investment to open hubs, the total transportation costs including customs, the external handling cost with seaport terminal or rail terminals and the internal handling cost within opened hubs (Eq 1):

, , , , , , , , , , , , , , , , , , , , , , , , | (Cos ) ( )
( ) _ int ( ) p m m t h h m o d o d h H m M t T o O d D m m t m t n m k t n m k t n n h o h h d h o h h d m M n N t T k K o O h H d D o h d Min t co Am z ct cd NV cm NV NV cm x x f ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ≠ ≠ = × × + + × + + × + + × + ×    Subject to: , , ( , , ) 1 m h z h m M h H h H z Dz Link h z m y z ∈ ∈ ∀ ∈ ∧ ∈ ∧ = =  1 2 2 1 , 2 1 2 ,

( , , ) 1

m h h h m M h H h H Link h h m y z ∈ ∈ ∀ ∈ ∧ = =  1 2 2 1 , 2 1 2 , ( , , ) 1 m h h h m M h H h H Link h h m y z ∈ ∈ ∀ ∈ ∧ = =  1 1 2 , , p {1, 2} (

, , ) 1 p h p h h h H P m

Link p h m y y z

∀ ∈ ∧ ∈ ∧ ∈ ∧ = + ≤ 1 2 1 2 1 2 1 2 1 2 2 3 , , 1

{2, 3} ,h Link( , h ,m)=1 Dist(h ,h ,m) D(m)

h h h h h m h H h y y z ∀ ∈ ∧ ∈ ∧ ∧ ≤ + ≤ 1 2 1 2 1 2 1 2 1 2 2 3 , , 2 {2, 3} ,h Link( , h ,m)=1 Dist(h ,h ,m) D(m) h h h h h m h H h y y z ∀ ∈ ∧ ∈ ∧ ∧ ≤ + ≤ 3 , Link( , ,3)=1 Dist(h,z,3) D(3) h z h h H z Dz h z y z ∀ ∈ ∧ ∈ ∧ ∧ ≤ ≤ 1 1, , , , , 1 p k k n t p p h h h H P k Kp t T x BigM z ∈ ∧ ∈ ∧ ∀ ∈ ∧ ∈ ≤ × 1 1 2 1 2 1 , , , , {2, 3 ( , , ) ( )
} p m k t p h h h t T m k K h H Dist h h m D m x BigM z ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∧ ≤ ≤ × 3, , , , ( , , 3) (3)
h p k t p h z h t T H z Dz k K Dist h z D x BigM z ∀ ∈ ∧ ∈ ∧ ∈ ∧ ≤ ∈ ∧ ≤ × , , , , , , m M n N t T k K o d ( , , ) 1 p k m k n t p m o d o d p P O D Link o d m x BigM y ∈ ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∈ ∧ ∈ ∧ = ≤ × ∧  , , , , , , m M t T k K o d ( , , ) 1 k p m m k n t p o d o d p P n N O D Link o d m y x ∈ ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∧ = ≤  , , , , , , , ( , ) , , , , ( , ) , , , , ' , , , ' h k

m k n t p m k n t p h p m k n t p h p p h h z h h m M m M z Dz h H

p P k K H t T x x x θ θ +∆ + + ∆ + ∈ ∈ ∈ ∈ ∈ ∀ ∈ ∧ ∈ ∧ ∈ = +   1 2 1 2 2 1 2 , , , ( , ) , , , , , , , , , h p m

k n t h h p m k n t p h h h z h H

z Dz m M m M k K t T H x x θ +∆ + ∈ ∈ ∈ ∈ ∀ ∈ ∧ ∈ ∧ ∈ =   (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11)
p m k t k t m h d z h d m M h H k K t T d D x ∈ ∈ ∀ ∈ ∧ ∈ ∧ ∈ = ×  , , , , , p m k t o d m t o d n k K m t T o O d D m M n N x NV CT f ∈ ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ≥ × (12)(13)(14)(15) , , , , , , D y



Constraints (2, 3, 4) express the single allocation of nodes to hubs. Constraints (5, 6, 7, 8) control the allocation mode to physical links. Outgoing hub flows exist only if the hub is active [START_REF] Rodriguez | Hub location under capacity constraints[END_REF][START_REF] Albareda-Sambola | Multiperiod Location-Routing with Decoupled Time Scales[END_REF][START_REF] Yang | Implications of postponement for the supply chain[END_REF] and require that this modal link should be already activated [START_REF] Van Hoek | The rediscovery of postponement a literature review and directions for research[END_REF][START_REF] Hsuan Mikkola | Supply-chain integration: implications for mass customization, modularization and postponement strategies[END_REF]. Constraints [START_REF] Schwartz | Designing distribution networks taking into account aspects of postponement, in The Development of Collaborative Production and Service Systems in Emergent Economies[END_REF][START_REF] Schulze | Location-allocation model for logistics networks with implementing commonality and postponement strategies[END_REF] ensure flow conservation at each period of time where Δ (o, h) +θ is the sum of the transports to h and transit time within h. Outgoing flows toward distribution centers must be equal to their respective demand [START_REF] Contreras | Stochastic uncapacitated hub location[END_REF]. Equation ( 17) computes the number of modal vehicles within the network where CTm is the capacity of a vehicle.

3.3.

Second Sub-problem: postponement location problem Given a set of located hubs Hl {h∈ H/žh=1} ∪ Dz and a set of active links L={(o, d, m); o∈ O, d∈ D, m∈ M/ ŷ m o,d=1}, we have to select the suitable location of postponement units in the designed distribution network. Location can be either on regional hubs, sub-regional ones or on local distribution centers, in order to minimize the total logistic costs (Equation 18) where b e h is a binary variable equal to 1 if the postponement unit is located on hub h at echelon e, while Bin (H, e) is a Boolean value equal to 1 if hub h is located at level e. ,

[ ( ) _ int ] [ ( ) + _ int b c e h h h Hl e E m m t m k t m o d m o d m o d m M t T o O d D m M t T k K o O d D m m t m o d m o d m m M t T o O d D m M t T k K o O d Min t postp co Am b ct cd cm NBV cm x ct cd cm NCV cm ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ = × × + + + × + × + + + × ×     , , , , , , ] m k n t p o d D p Pk x ∈  Subject to: h ( , ) 1 e h h e E b z H Bin H e ∈ = ∀ ∈ ∧ =  ) 1 3 1 2 3 , , 3 , , 1 1 , 1 , , , , , , , , , , , , 2 2 , 2 , , , , , , , , , (Cos _ )
x BigM b ,

k m k t m h h h h h b m M p P y t T k K h h Hl ∈ ∈ ≤ × × ∀ ∈ ∧ ∈ ∧ ∈  1 2 1 2 2 3 , , 2 3 , , 1 2 3 , x BigM (b + b ) t T k K , , k m k t h h h h h h b m M p P y h h h Hl ∈ ∈ ≤ × × ∀ ∈ ∧ ∈ ∧ ∈  2 2 2 2
, , ', ( , ) , , ,,,' ,' , ,

( , ') / , , , c k m n k t h z m k t m h z h z p h z z m M m M k K p P b x y x Conv n Cs k k b t T h Hl z Dz k K θ +∆ + ∈ ∈ ∈ ∈ = × × × × ∀ ∈ ∈ ∈ ∈   ) 1 2 1 2 1 2 2 2 1 ', , ', ( , ) , , ' 1 , , , 3 , , , , , ' , 2 ( )
1 2 ( ) ( , ') (1 )

/ , , , c k m n k t h h m k t m m h h h h h z p h z

h m M m m M k K p P z Dz b x y y x Conv n Cs k k b t T h h Hl k K θ +∆ + ∈ ∈ ∈ ∈ ∈ = × × × × × - ∀ ∈ ∈ ∈   ) ) (16) (17) (18) (19) (20) (21) (22) (23) 1 2 1 1 1 1 , , ', ( , ) , , 1 , , , , , ' , , 1 2 ( ) ( , ') (1 )
/ , , ,

c k m n k t h h m k t m h z h z p h z h m M m M k K p P z Dz b x y x Conv n Cs k k b t T h h Hl k K θ +∆ + ∈ ∈ ∈ ∈ ∈ = × × × × - ∀ ∈ ∈ ∈   ) , , , , , , 3 1 ,
; b k m k t p o d m t o d m k K p P m x NBV NV t T m M h Hl h Dz CT ∈ ∈ ≥ × ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈  3 3 , , , , , , , 3 2 , , [ ] / c
k m k n t p h h m t e h h h n e E k K p P n N m x NCV b t T m M h Hl h Dz CT f ∈ ∈ ∈ ∈ ≥ × ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ×   2 3 2 3 1 2 , , , , , , , 1 2 , 1 2 3 2 , ( ) [ ]
/ , Dz c k m n k t p h h m t h h h h n k K m p P n Pg x NCV b b t T m M h h Hl h CT f ∈ ∈ ∈ ≥ + × ∀ ∈ ∧ ∈ ∧ ∈ ∧ ∈ ×  1 1 1 2 1 2 2 1 2 , , , , 1 , , ', " , , , , , ' 1 1 2 2 , , , 2 ( ) [( ( , , ) ( )) ( , , ) ( ', , ) ( ) ( )) ]
[ ( ", , )

k t n k t n m z h postp z h z n m m m M h h Hl k t n k t m z h h postp z h z h n k t n k z h postp z n z Dz D z m h z D Conv n y f D m h h m h z D Conv n y b f D m h z D f δ θ δ δ θ δ θ ∈ ∈ ∀ ∈ ∆ = ∆ + × + × × × + ∆ + ∆ + + × + × × × × ∆ + × + ×  2 2 , 2 " , , 3 ( )]
( )

t m h h z k t postp z z z Conv n b y D Conv n b L θ × × × + × × × ≤ ∆
Constraint [START_REF] Gelareh | Liner shipping hub network design in a competitive environment[END_REF] translates that a postponement activity can only be located on activated hubs, at only one level. "Continuous flows" x m,k,t o,h exist only if postponement units are located before that hub if a modal link is activated (20, 21). (22, 23, 24) express a flow balance at each period and on each and computes the ingoing discrete flows to each hub depending on location of postponement units. Conv(n) is the conversion ratio from unit of packaging product under n commodity to a continuous unit (tons for example) and Cs(k, k') is the amount of component k needed to produce a unit of product k'. (25, 26, 27) assess the number of bulk and container vehicles within the network. Let ΔLz be the requirement of a market zone on the level of service. (28) expresses the requirement on the service level where Δ(m', h, z) is the transportation time from postponement units, h δ the transit processing time within transited hubs and θpostp the unit postponement time.

Illustrative study

This case study aims to illustrate the application of the proposed models. It concerns the location of blending units within East Africa for specific industries involving hybrid (discrete-continuous) flows, like the fertilizer industry. Three production zones located in Morocco, Ethiopia and Nigeria and considered. Five Regional hubs are defined: Kenya-Angola-Tanzania-Djibouti, so that ten sub-regional hubs: Nairobi-

(25)

Kisumu-Dodoma-Arusha-Tabora-Kuito-Tete-Lichinga-Kigali. The considered data are summarized in Tables 1 and2. The demand (aggregated on one year) varies for each zone. We assume that shipment is done every month and that all the market zones require the same service level. Distances and traveling times are extracted from Google Maps. The problem was solved using the Xpress-IVE solver tools. The results of models 1 and 2 are summarized in Tables 3 and4.

Conclusion et research perspectives

In the context of mass customization, postponement activities may provide an answer to fulfill customized orders, increase customer responsiveness and increase service level. However, the literature combining design of logistic hub networks and implementation of postponement facilities is scarce and usually assumes that the location of the distribution centers are already known. In order to address the problem, we have developed and tested a two-phase deterministic mathematical programming model where, as a first step, we design incapacitated discrete logistic hub, then allocate postponement services on some hubs. In our future work, this model will be coupled with a discrete event simulation model in order to take into consideration uncertainties on the demand and on the availability of the resources. Simulation will also allow to assess postponement capacities and to refine logistic costs, these results being re-injected in the mathematical model as new constraints.

Table 1 .

 1 List of transport cost parameter

Table 2 .

 2 List of other cost parameters

	Parameter	Discrete flow	Continuous flow	Cost		Regional hub	Sub-regional hub
	Transport capacity-	50 palette	30000 Tons	Hub Location cost	5760000	2880000
	sea			($/year)
	Transport capacity-	40 palette	25000 Tons	Blending location 30000	30000
	rail			Intern	handling	15$/Tons	35$/Tons
	Transport capacity-	40 palette	25000 Tons	(bulk)
	road			Extern	handling	100 $/train	120 $/train
	Rail unit transport	0.06/train/km	0.04/Train/km	(rail)
	cost			Intern	handling	5 $/palette	10 $/palette
	Road unit transport	3.75$/truck/km 2$/truck/km	(discrete)
	cost			Extern	handling	40 $/train	50 $/train
				(rail) (discrete)

Table 3 .

 3 Results of sub-model 1

	Located hub	Plant	Sub-hub	% DC	Total logistic cost 216 M$
	Kenya	Morocco	Kisumu-	25%	Sea Transport	9.52 M$
		Ethiopia	Lichingua		Rail Transport	8.75 M$
	Mozambique	Morocco	Tete	10%	Road Transport	21.97 M$
	Angola	Nigeria	luau	10%	Intern handling	109.45 M$
	Tanzania	Morocco	Kuito-Kigali	10%	Extern handling	26890.5 $
		Ethiopia			Location	66.3 M$

Table 4 .

 4 Results of sub-model 2

	Blending	Kenya-Mozambique-Angola-	Total logistic cost	157 M$
	location	Tanzania	Total bulk costs	47.1 M$
	Level	1	Total discrete costs	109.9 M$