domain: sde sdv

text:
Maryem Hamidi
email: hamidi.maryem@gmail.com

Tarik Bouramtane

Shiny Abraham

Ilias Kacimi

Laurent Barbiero

Nadia Kassou

Vincent Valles

Gad Levy

Landslide hazard assessment in the heterogeneous geomorphological and environmental context of the rif region, morocco -A machine learning approach

Keywords: landslide, machine learning, Morocco, susceptibility mapping, remote sensing, Rif

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

INTRODUCTION

Landslides pose a significant danger in various mountainous regions around the globe and can result in devastating consequences such as loss of property, economic instability, and loss of life [START_REF] Prasad | Natural Disasters[END_REF]. The destructive impact of landslides is increasing globally due to deforestation, global population growth, and climate change, whereby the latter alters the stability of slopes and bedrock through changes in precipitation and/or temperature [START_REF] Saha | Hybrid ensemble machine learning approaches for landslide susceptibility mapping using different sampling ratios at East Sikkim Himalayan, India[END_REF]. Landslides can be local and superficial but could also occur on a large scale and suddenly, affecting the slopes of entire catchments and leading to severe landscape damage. Therefore, understanding the dynamics of landslides and developing effective mitigation measures are essential for minimizing the impact of these hazard events and ensuring the safety of vulnerable populations and infrastructure.

Africa, notably, stands out as a significant landslide-prone area due to high population growth and the influence of climate change, resulting in increased precipitation and more frequent and violent rainstorms [START_REF] Broeckx | A data-based landslide susceptibility map of Africa[END_REF][START_REF] Gariano | Landslides in a changing climate[END_REF]. The tectonically active region extending from Morocco to Egypt contributes to the heightened susceptibility of the northern margin of Africa to landslides [START_REF] Poggi | A probabilistic seismic hazard model for North Africa[END_REF]. Despite the critical necessity for landslide susceptibility mapping in North Africa, only a few studies have been conducted, and the majority of them have focused on Morocco [START_REF] Abdı | Landslide Susceptibility Mapping Using GIS-based Fuzzy Logic and the Analytical Hierarchical Processes Approach: A Case Study in Constantine (North-East Algeria)[END_REF][START_REF] Anis | GIS-based landslide susceptibility mapping using bivariate statistical methods in North-western Tunisia[END_REF]. Integrating landslide susceptibility mapping into the disaster management plans of each exposed country is crucial for improving disaster management and mitigation strategies (P. Zhao et al., 2022).

The Taounate province in Northern Morocco is prone to numerous landslides and is located within the geologically diverse Intra and Mesorif domains, as well as the geologically consistent Prerif region. A recent study has identified over 1,000 landslide cases in Taounate alone, equivalent to 25 movements/km2 [START_REF] Abidi | Landslide Susceptibility Assessment Using Evidence Belief Function and Frequency Ratio Models in Taounate city (North of Morocco)[END_REF]. The geology of this area, characterized by active tectonics and rugged terrain, and its climate, characterized by long periods of drought followed by periods of heavy precipitation in the form of extreme events, are considered the leading causes of this natural hazard [START_REF] Benabdelouahab | Analysis and trends of rainfall amounts and extreme events in the Western Mediterranean region[END_REF]. Moreover, all the studies conducted to generate landslide susceptibility maps are case-by-case studies and do not cover the entire province. For instance, investigations by Jemmah & Brahim, 2018, and the research conducted by [START_REF] Benchelha | Landslide Susceptibility Mapping in the Commune of Oudka, Taounate Province, North Morocco: A Comparative Analysis of Logistic Regression, Multivariate Adaptive Regression Spline, and Artificial Neural Network Models[END_REF], focused on the province's northern regions, where they performed a comparative assessment of various techniques for mapping landslide susceptibility. Nevertheless, previous research efforts have not adequately addressed the overall complexities of the entire region under study. Producing clear maps is essential in mapping landslide susceptibility, and various techniques have been suggested and utilized for this aim, including qualitative and quantitative approaches. Qualitative methods require an expert to estimate landslide potential and slope movement. On the other hand, quantitative methods, such as statistical or probabilistic techniques, have proven effective in identifying areas that require intervention or development and are used by policymakers with input from scientists, engineers, and the general population (Bravo-López et al., 2022; [START_REF] Reichenbach | A review of statisticallybased landslide susceptibility models[END_REF].

Numerous statistical techniques have been used in order to generate landslide susceptibility maps in the Moroccan Rif region. For example, Hierarchical Fuzzy Inference Systems (HFIS) that were applied in the Rif Mountains' center region [START_REF] Ozer | On the use of hierarchical fuzzy inference systems (HFIS) in expert-based landslide susceptibility mapping: The central part of the Rif Mountains (Morocco)[END_REF] and fuzzy Analytical Hierarchy Process (FAHP) were combined for landslide susceptibility modeling in southeast Morocco (Sadiki et al., 2023); Weight Of Evidence (WOE) was applied in Tetouan-Ras Mazar [START_REF] Elmoulat | Landslides susceptibility mapping using GIS and weights of evidence model in Tetouan-Ras-Mazari area (Northern Morocco)[END_REF], an area of the Rif chain (Es-smairi et al., 2021) and Taounate-Ain Aicha (Jemmah & Brahim, 2018). Logistic regression (LR) was used in Taounate-Oudka and the Sahla catchment area (El-Fengour et al., 2021), while the Analytical Hierarchical Processes (AHP) method was applied in Tangiers [START_REF] Brahim | Landslide susceptibility mapping using AHP method and GIS in the peninsula of Tangier (Rif-northern morocco)[END_REF], the Larache Province [START_REF] Hamdouni | The Prevention of Landslides Using the Analytic Hierarchy Process (AHP) in a Geographic Information System (GIS) Environment in the Province of Larache, Morocco[END_REF] and parts of the Rif mountains (Es-smairi et al., 2021). However, these methods rely on assumptions about the data, such as normality and linearity, which can limit their accuracy in analyzing complex interrelationships between different causative factors. Furthermore, reclassifying continuous causative factors to prepare the dataset can be challenging and affect the nature of the data [START_REF] Chen | Landslide Susceptibility Modeling Based on GIS and Novel Bagging-Based Kernel Logistic Regression[END_REF], resulting in decreased model accuracy. Additionally, classical statistical methods, which require large sample sizes and are less able to process high-dimensional and complex data (He et al., 2012), The geo-environmental processes controlling landslides in the Taounate province are complicated and remain unknown, highlighting the need for further research and effort in landslide susceptibility modeling and mapping. However, there are opportunities for improvement by leveraging the strengths of high-resolution, freely available remote sensing data and the efficiency of advanced machine learning methods. Machine learning techniques have lately gained significant popularity in assessing landslide susceptibility [START_REF] Rahmati | Land subsidence modelling using tree-based machine learning algorithms[END_REF]. The accessibility of satellite imagery has enhanced machine learning's success, geographic information systems (GIS), and the scarcity of historical landslide records [START_REF] Ado | Landslide Susceptibility Mapping Using Machine Learning: A Literature Survey[END_REF]. Despite the numerous advantages offered by machine learning, including its ability to handle high-dimensional data efficiently, fit nonlinear relationships between targets and factors (X. [START_REF] Zhou | Geoinformation-based landslide susceptibility mapping in subtropical area[END_REF], as well as its adaptability, reproducibility, and ability to quantitative analysis of factors contributing to landslide development. The use of machine learning methods for landslide susceptibility mapping has been limited in North Africa, particularly in Morocco.

Few research studies have evaluated the efficiency of Artificial Neural Networks (ANN) and LR in identifying suitable methods for limited regions within the Taounate province. In a study by [START_REF] Benchelha | Landslide Susceptibility Mapping in the Commune of Oudka, Taounate Province, North Morocco: A Comparative Analysis of Logistic Regression, Multivariate Adaptive Regression Spline, and Artificial Neural Network Models[END_REF], the Multivariate Adaptive Regression Spline (MARSpline) model demonstrated a higher success rate compared to LR and ANN models. However, it should be mentioned that this study was verified using a single assessment measure (Area Under Curve (AUC)) and was based on eight conditioning factors. In other research, Sahrane et al., 2023 used LR and ANN models to examine the effect of landslide inventory maps on landslide susceptibility in two distinct geomorphologically different regions in Taounat province. In the Atlas Mountains, the Support Vector Machine (SVM) model outperformed the Weight of Evidence (WoE) and Radial Basis Function Network (RBFN) models despite being utilized only once in Morocco [START_REF] Naceur | Performance assessment of the landslide susceptibility modelling using the support vector machine, radial basis function network, and weight of evidence models in the N'fis river basin, Morocco[END_REF]. The case studies mentioned have their limitations, as the effectiveness of machine learning techniques can be influenced by the unique geo-environmental characteristics of each region, various contributing factors, and the accuracy of input data (P. Zhao et al., 2022). Hence, to ensure the generation of precise landslide susceptibility maps, it is essential to assess and contrast alternative advanced machine learning approaches capable of adequately addressing the overall study area complexities when performing landslide susceptibility mapping.

Addressing the current limitations of models used and the data gaps in Morocco's Rif mountains, we concentrated on analyzing the performance of four machine learning methods: SVM, LR, Multivariate Discriminant Analysis (MDA), and Classification and Regression Trees (CART) for landslide susceptibility mapping. These models were chosen based on a literature review and the recognition that different linear and nonlinear factors influence landslide occurrence. The aim of this research study is not limited to the assessment of machine learning models but also to develop an extensive understanding of the geo-environmental factors controlling landslides in one of the most affected and complex regions of North Africa. Moreover, the strength of this study lies in its extensive and significant database, comprising 255 landslides and 15 geo-environmental factors using remote sensing and fieldwork data. These contributions represent a substantial addition to the scientific community.

STUDY AREA

The geographical area studied covers a mountainous province of around 5616 km² in the Fez-Meknes region of North Morocco (Figure 1). In 2009, the province's population was approximately 678,000, giving it an average density of 121 inhabitants per km², relatively high for a mountainous area. The province straddles three structural subdomains of the Rif belt: Intra-Rif, Meso-Rif, and Pre-Rif, and therefore is characterized by different tectonic units and a high diversity of lithological units (limestone, marl, flysch, etc.) (

DATA AND METHODS

A schematic representation of the application procedure is shown in Figure 2. The main steps in the methodological approach comprise data collection and organization, the model's development and validation, and the contribution of conditioning factors.

Data collection

Generating landslide inventory map

Past and present landslides are critical indicators of future landslides (Guzzetti et al., 1999). Therefore, an inventory of landslides, based on historical data collected through exhaustive investigations by Maurer, 1968 and high-resolution Google Earth images, was utilized to create a map of landslides in the study region. Google Earth Pro software was employed to validate the existence of historical data and collect new landslide and non-landslide points. Non-landslide locations were randomly selected from areas without landslide, ensuring a minimum distance of 2000 meters from landslide sites. In order to produce a single landslide inventory map, six different types of landslides, including rockfalls, debris flows, and complicated landslides, were identified. The inventory map contained 255 landslides and 255 non-landslide points in Taounate province (Figure 1). Figure 2 illustrates the methodological process used in this study.

Generating landslide conditioning factors

The interaction of different factors causes the occurrence of a landslide. There are no specific guidelines for selecting these factors, each of which may contribute to a greater or lesser extent in different regions. However, it is commonly assumed that landslides will happen under similar conditions as in the past. [START_REF] Youssef | Landslide susceptibility mapping using machine learning algorithms and comparison of their performance at Abha Basin, Asir Region, Saudi Arabia[END_REF]. Based on the literature assessment and the study area's environmental and geographical characteristics, fifteen landslide conditioning factors were chosen. They can be divided into four categories (Table 1, Figure 3), namely topographic (altitude, Topographic Position Index (TPI), Topographic Wetness Index (TWI), slope, aspect, curvature, profile curvature, plan curvature), hydrological (precipitation, distance to river), geological (lithology, distance to faults), and land cover related (distance to roads, Land Cover and land Use (LCLU) and Normalized Difference Vegetation Index (NDVI)). The discrete variables were distributed following the properties of the data, while the topographic factors were derived using a 30 m resolution Digital Elevation Model (DEM). Using the GIS tool, all conditioning factors were converted into a raster format and synchronized at a spatial resolution of 10×10 m cell size.

Slope aspect

The slope's aspect or exposure indicates the direction in which it is oriented [START_REF] Barman | Application of index of entropy and Geospatial techniques for landslide prediction in Lunglei district, Mizoram, India[END_REF]. Conditioning exposure to sunlight and drying winds that control soil moisture and evapotranspiration [START_REF] Tang | Evaluating landslide susceptibility based on cluster analysis, probabilistic methods, and artificial neural networks[END_REF]) is an indicative factor of landslide susceptibility. This factor can influence other factors, including wind direction, rainfall regime, amount of sunlight, hydrology, vegetation development, and soil moisture [START_REF] Xing | Large-Scale Landslide Susceptibility Mapping Using an Integrated Machine Learning Model: A Case Study in the Lvliang Mountains of China[END_REF].

Curvature

The curvature is the change in the topographic surface's slope over a narrow curve arc, also called slope shape or topographic curvature. It influences the direction of surface flow and the constraint on the slope material.

Pro le curvature

Profile curvature is the term used to describe the curvature in the primary direction of the steepest slope. This factor impacts how quickly and slowly the flow moves, impacting erosion and deposition [START_REF] Xie | A Novel Hybrid Method for Landslide Susceptibility Mapping-Based Geo-Detector and Machine Learning Cluster: A Case of Xiaojin County, China[END_REF]. A positive or negative value denotes the cell is part of a concave or convex upward slope.

Plan curvature

The curvature of a plane is defined by the contour line created when the surface intersects a horizontal plane. It controls the water direction during downslope flow, which can condition landslide events (Habumugisha et al., 2022; Tavakkoli Piralilou et al., 2019). The cell is a portion of a laterally concave or convex slope, whether the value is positive or negative.

TPI

The Topographic Position Index (TPI) is a parameter to describe the terrain (Andrew D. Weiss, 2001). It refers to the difference in altitude between a point and its surroundings. It is the outcome of comparing each cell's height in a DEM to the average elevation of its surrounding cells. The TPI is positive for cells elevated above their surroundings (ridges and hilltops) and negative for depressions.

TWI

The water saturation zones and soil moisture content are represented by the Topographic Wetness Index (TWI). This factor was chosen because it explains how topography affects soil moisture. Higher values of TWI can be related to higher risks of landslides [START_REF] Meena | Assessing the importance of conditioning factor selection in landslide susceptibility for the province of Belluno (region of Veneto, northeastern Italy)[END_REF].

Lithology

Lithology is one of the most frequently employed parameters in landslide susceptibility studies due to its impact on the geo-mechanical properties of the terrain (Costanzo et al., 2012). The variation of lithological units often leads to differences in permeability, density, hardness, and strength of soils and rocks. The lithological map was created from various sources, including the United States Geological Survey's (USGS) geological data for Africa and the 1:200,000 geological map of Morocco. Eight lithological units were distinguished and assigned the following codes:

• Quaternary (marls, sandstone, limestone, and alluvium) = 0 • Lower Jurassic (limestone and molasses), Cretaceous (marls and marly limestone), and Cretaceous-Tertiary (marls, marly limestone, and flysch) = 1 • Tertiary (limestone and flysch) =2 • Lower Cretaceous (Flyschs), tertiary intrusive rock, tertiary intrusive rock, and Triassic (basalt) = 3.

Precipitation

Precipitation is one of the main initiators or triggers of landslide risk. It raises pore pressure and water content of the soil, reducing soil cohesion and potentially leading to slope instability and mass movement [START_REF] Habumugisha | Landslide Susceptibility Mapping with Deep Learning Algorithms[END_REF]. The calculation of monthly average precipitation spanning 16 years (2005-2021) was performed using data from the Climate Hazards Group Infra-Red Precipitation with Stations (CHIRPS) [START_REF] Funk | The climate hazards infrared precipitation with stations-A new environmental record for monitoring extremes[END_REF]. The CHIRPS data were preferred because several studies concluded that CHIRPS performs better than TRMM 3B42 (Tropical Rainfall Measuring Mission version 7) and PERSIANN-CCS (Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks -Cloud Classification System) [START_REF] Karmouda | Geo-statistical and hydrological assessment of three satellite precipitation products over Ouergha basin (Northern Morocco)[END_REF]. For the past 16 years, the research region has received an average of 35.19 to 68.50 mm of rain each month.

Land cover/ land use

The physical characteristics of soils, as well as the hydrological processes of rainfall distribution, infiltration, slopes, and runoff, may all be influenced by land use [START_REF] Meena | Assessing the importance of conditioning factor selection in landslide susceptibility for the province of Belluno (region of Veneto, northeastern Italy)[END_REF]. These data for 2021 in our study area were directly extracted from ESRI's annual 10m land cover maps, derived from ESA's Sentinel-2 imagery at 10m resolution. These were produced by the Impact Observatory, which employed a deep-learning model trained for land classification using billions of human-labeled pixels. The methodology for encoding the land cover parameters was as follows: Water = 1, Built area = 2, Bare ground = 3, Range Land = 4, Flooded vegetation = 5, Crops = 6, and Trees = 7.

Distance to roads

One of the human factors that could promote various forms of erosion is the distance to roads, especially the occurrence of landslides. Road construction requires excavation and dumping of material, which impacts slope stability, and vibrations generated by human activities can trigger mass movements [START_REF] Vakhshoori | Landslide susceptibility mapping by comparing weight of evidence, fuzzy logic, and frequency ratio methods[END_REF]. The road distance calculation was derived from the road network extracted from Google Earth through a GIS tool.

NDVI

The presence and type of vegetation cover are essential to mass movement, influencing slope stability by water absorption and reinforcing soil cover [START_REF] Xing | Large-Scale Landslide Susceptibility Mapping Using an Integrated Machine Learning Model: A Case Study in the Lvliang Mountains of China[END_REF]. A high (NDVI) value may also suggest an elevated likelihood of such events [START_REF] Huang | Landslide Susceptibility Prediction Considering Regional Soil Erosion Based on Machine-Learning Models[END_REF]. Vegetation density was obtained by calculating the Normalized Difference Vegetation Index (NDVI) on Google Earth Engine (GEE), using a one-year time series (January to December 2021) derived from cloudfree images of the Sentinel-2A and Sentinel-2B satellites. The NDVI map was calculated according to Equation 1:

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁 (1)
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 (4) (1)
NIR and RED refer to the spectral reflectance measurements obtained within the red and nearinfrared bands of the electromagnetic spectrum, respectively. The NDVI value varies within the range of -1 to 1.

Distance to faults

The distance from a fault is a significant geological factor that can potentially trigger landslides, and it has been utilized in numerous prior studies related to landslide susceptibility (Habumugisha et al., 2022). Areas near geological faults are at high risk for landslides because the surrounding rock experiences significant fracturing, decreasing its strength due to tectonic fractures [START_REF] Chen | Landslide spatial modeling: Introducing new ensembles of ANN, MaxEnt, and SVM machine learning techniques[END_REF].

Distance to rivers

Many landslides are caused by rivers because of their power to erode exposed material, thus becoming a key factor in landslide susceptibility (Pradhan et al., 2010). This factor was computed using the "Euclidean Distance" from the drainage network derived from the SRTM DEM data and refined using Google Earth for improved accuracy.

Modeling using Machine Learning Algorithms

Multivariate discriminant analysis

Multivariate discriminant analysis (MDA) is a multivariate statistical method that complements linear discriminant analysis (LDA). It performs a mixture of classifiers for clustering by merging the implementation of simple linear and non-linear combination models. The MDA algorithm has many advantages, including its ability to handle large and high-dimensional datasets, avoid overfitting of datasets, and be easy to use (Hastie & Tibshirani, 1996). It is a widely employed supervised machine learning methodology utilized for both regression and classification tasks in various fields, including natural resources modeling and landslide susceptibility assessments [START_REF] Kalantar | Optimized Conditioning Factors Using Machine Learning Techniques for Groundwater Potential Mapping[END_REF].

Logistic regression

Logistic regression (LR) is a simple, rapid, and highly efficient classification method (Rajaneesh et al., 2022). It is one of the most commonly utilized machine learning algorithms, originally developed within the statistical field (Cabrera, 1994). The generalization of the linear model relies on employing the logistic function to represent a binary response variable (dependent variable, landslide occurrence) and numerous explanatory variables (independent variables, conditioning factors) to predict the likelihood of specific events, such as landslides [START_REF] Xie | A Novel Hybrid Method for Landslide Susceptibility Mapping-Based Geo-Detector and Machine Learning Cluster: A Case of Xiaojin County, China[END_REF].

The classification output is between 0 and 1, making it applicable for probabilistic interpretation in several study areas [START_REF] Qu | Review of bankruptcy prediction using machine learning and deep learning techniques[END_REF].

Classi cation and regression Ttees

Classification and regression trees (CART), introduced by Breiman in 1984, represent predictive algorithms employed in the machine learning field. Instead of using stopping rules, CART grows an extensive tree to produce a sequence of subtrees for solving classification problems (Loh, 2011). The decision tree algorithm has been demonstrated to be highly effective in dealing with classification. It has the noticeable advantage of handling small datasets and scaling to significant problems (Markham et al., 2000). The model is a non-parametric procedure for predicting output variables using input variables [START_REF] Prakash | Machine Learning Methods of Kernel Logistic Regression and Classification and Regression Trees for Landslide Susceptibility Assessment at Part of Himalayan Area, India[END_REF]. CART has found extensive application in landslide studies across various global regions (Felicísimo et al., 2013; Rabby & Li, 2020; Vorpahl et al., 2012), although, to the best of our knowledge, it has not been applied in Morocco. In the context of landslide prediction, the CART algorithm undergoes a four-step training process: (i) tree construction, (ii) determination of tree construction termination, (iii) tree pruning, and (iv) selection of the optimal tree for classifying landslide and non-landslide categories (Loh, 2011).

Support Vector Machine

Support Vector Machine (SVM) is an efficient and robust supervised technique applicable for both regression and classification tasks. (Cortes & Vapnik, 1995; Cristianini & Shawe-Taylor, 2000). The most crucial feature of SVM is its kernel function, which transforms the input variables into a complex and multi-dimensional linear space through nonlinear transformation. This generates an optimal separation hyperplane between two classes [START_REF] Qu | Review of bankruptcy prediction using machine learning and deep learning techniques[END_REF] to maximize the margin between them, thereby making it a nonlinear classifier. There are typically four categories of Kernel functions utilized in SVM: linear, sigmoid, radial basis, and polynomial functions. In this study, the linear function was chosen for the SVM model to compute the landslide susceptibility index. SVM classification has demonstrated its effectiveness in landslide susceptibility mapping, outperforming various other machine learning methods [START_REF] Xie | A Novel Hybrid Method for Landslide Susceptibility Mapping-Based Geo-Detector and Machine Learning Cluster: A Case of Xiaojin County, China[END_REF], especially in mountainous areas.

Model validation and performance

As landslide mapping involves binary classification (absence or presence of landslides), it is crucial to assess the quality of landslide probability estimates. This assessment should consider prediction accuracy, misclassification errors, and outcomes [START_REF] Merghadi | Machine learning methods for landslide susceptibility studies: A comparative overview of algorithm performance[END_REF]. Before model validation, the dataset must be divided into training and testing parts. Previous studies have commonly used simple random division [START_REF] Kadavi | Application of Ensemble-Based Machine Learning Models to Landslide Susceptibility Mapping[END_REF]. Nevertheless, we chose to employ robust statistical validation through k-fold cross-validation to address the issue of imbalanced landslide samples and mitigate bias [START_REF] Pal | Data Classification with k-fold Cross Validation and Holdout Accuracy Estimation Methods with 5 Different Machine Learning Techniques[END_REF]. Using this approach, the dataset is divided randomly into 10 subsets or folds, with each fold employed to validate the model trained on the remaining folds. This holdout process is repeated 10 times. The evaluation is intended to identify and compare the accuracy of the four machine learning models across the Taounate province. The landslide susceptibility values of the four models, which ranged from 0 to 1, were compared with the landslide feature based on two predicted landslide classes: "positive" if the value was more significant than 0.5 and "negative" if the value was less than 0.5. From these values, six statistical measures were calculated, including precision, sensitivity (recall), specificity, accuracy, F1-measure (or F1-score), and receiver operating characteristics (ROC), as follows:

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁 (1)
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (4) 𝑇𝑇𝐴𝐴𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑃𝑃𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (5) 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (6) DACCi = 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 -𝐴𝐴𝐴𝐴𝐴𝐴 𝑖𝑖 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 × 100 (7) (2)

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁

(1)

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (4) 𝑇𝑇𝐴𝐴𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑃𝑃𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (5) 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (6) DACCi = 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 -𝐴𝐴𝐴𝐴𝐴𝐴 𝑖𝑖 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 × 100 (7) (3)

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁

(1)

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (4) 𝑇𝑇𝐴𝐴𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑃𝑃𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (5) 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (6) DACCi = 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 -𝐴𝐴𝐴𝐴𝐴𝐴 𝑖𝑖 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 × 100 (7) (4)
𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁 (1)
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (4) 𝑇𝑇𝐴𝐴𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑃𝑃𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (5) 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (6) DACCi = 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 -𝐴𝐴𝐴𝐴𝐴𝐴 𝑖𝑖 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 × 100 (7)
(5)

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝑁𝑁𝑁𝑁𝑁𝑁 -𝑁𝑁𝑅𝑅𝑁𝑁 𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑁𝑁𝑅𝑅𝑁𝑁 (1)
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 (2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (3) 𝐴𝐴𝑆𝑆𝑆𝑆𝐴𝐴𝐴𝐴𝐴𝐴𝑆𝑆𝑆𝑆 = 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (4) 𝑇𝑇𝐴𝐴𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑃𝑃𝑆𝑆 = 𝑇𝑇𝑇𝑇 𝑇𝑇𝑁𝑁 + 𝐹𝐹𝑇𝑇 (5) 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 + 𝐹𝐹𝑇𝑇 (6) DACCi = 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 -𝐴𝐴𝐴𝐴𝐴𝐴 𝑖𝑖 𝐴𝐴𝐴𝐴𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 × 100 (7) (6)
where: "P" represents the total number of landslides, "N" represents the number of nonlandslide points, "True Positives (TP)" and "True Negatives (TN)" denote correctly classified samples, and "False Positives (FP)" and "False Negatives (FN)" represent misclassified samples.

The ROC curve, as introduced by Bradley, 1997, is a commonly used method for assessing the predictive capabilities of a model. The area under the ROC curve (AUC) serves as a performance metric that considers specificity on the x-axis and sensitivity on the y-axis [START_REF] Jiang | Leveraging Machine Learning and Simulation to Advance Disaster Preparedness Assessments through FEMA National Household Survey Data[END_REF]. AUC values are between 0 and 1 and can be classified as excellent (0.9-1), very good (0.8-0.9), good (0.7-0.8), average (0.6-0.7) and poor (0.5-0.6). According to Mathew et al., 2009, the AUC is a noteworthy metric and one of the most valuable precision statistics for landslide susceptibility analysis.

Factor importance

For geo-environmental models, it is necessary to determine which factors are essential or contribute the most to the disaster. There is no universal guide for identifying these factors, and several methods have been proposed in recent years. We opted for a jack-knife test based on removing one factor at a time and using the remaining factors for modeling ( ACC all represents the calculated model accuracy when all parameters are considered. ACC i denotes the accuracy value of the model when an indicator 'i' is excluded from the input data, and DACC i represents the corresponding percentage decrease in accuracy.

RESULTS

Comparative performance analysis of machine learning models

According to the AUC index, the variation in model performance was minimal. All four models performed excellently for landslide predictive ability with AUC values above 0.95. The CART model achieved the highest score (AUC=0.971), followed by LR (AUC=0.963), with SVM (AUC=0.961) and MDA (0.954) following suit (Figure 4). Evaluating model performance relying only on one assessment measure, such as AUC, might not be suitable, as a high AUC value does not invariably signify elevated accuracy in spatial prediction [START_REF] Kalantar | Optimized Conditioning Factors Using Machine Learning Techniques for Groundwater Potential Mapping[END_REF]. Therefore, five statistical measures were used: precision, sensitivity, specificity, accuracy, and F1-measure. All the results are presented in Table 2. The SVM model exhibited superior performance in classifying both landslide and nonlandslide in both the training and test datasets,

Landslide susceptibility maps

Figure 5 illustrates the landslide susceptibility maps. The use of the four models resulted in a consistent distribution of susceptibility prediction categories, with high susceptibility areas primarily concentrated in the northern and northwestern regions of the study area. Conversely, regions characterized by low susceptibility were situated in the southern part. Meanwhile, areas displaying varying degrees of susceptibility, ranging from high to moderate, were observed in the central, east-central, and west-central regions of the study area. Subsequently, the relative proportions of the four classes in each model were computed (Figure 6). The proportion of the very high susceptible area was comparable between the CART and SVM models, constituting 42% in both cases, slightly higher for the LR model (45%) and lower for the MDA model (40%). On the other hand, the differences between the models were stronger for the low susceptibility category, ranging from 28% for the CART model to 37% for the MDA model. For the medium and high susceptibility categories, there were some similarities between the four models, as they occupy less than 16%.

Contribution of conditioning factors to landslide distribution

Figure 7 illustrates the explanatory variables' contribution for each model and is represented by the number of observations appearing on the x-axis. Precipitation had very high DAAC values ranging from 6.52% to 14.66%. From this finding, it is concluded that precipitation is the most crucial factor that significantly affects the distribution of landslides for each classifier in the research region, followed by distance to the rivers with DAAC values ranging from 4.21% to 3.04%. However, it appears that distance to roads and land cover/land use held relatively lower significance within the study area, as indicated by the outcomes of the four classifiers. The comparative importance of precipitation, distance to rivers, land cover, and distance to roads was similar among the four models. However, the relative importance of the remaining conditioning factors varied significantly across the models. Notably, the SVM model attributed the highest DAAC value to precipitation (14.66%), marking it as the most influential factor. It was followed by slope (5.25%), distance to rivers (4.16%), altitude (2.84%), and distance to faults (2.84%).

DISCUSSION

Factors contributing to landslide risk in the Rif region

The investigation into landslide susceptibility is essential in land management, planning, and the development of hilly regions. These susceptibility maps are valuable tools for government agencies, offering essential guidance for effective decision-making in landslide hazard and risk analysis, as well as for general planning and assessment of landslides [START_REF] Roccati | Heavy Rainfall Triggering Shallow Landslides: A Susceptibility Assessment by a GIS-Approach in a Ligurian Apennine Catchment (Italy)[END_REF]. This study is innovative as it marks the initial endeavor to produce a comprehensive landslide susceptibility map encompassing the entire study region. The results revealed that nearly half, approximately 40% to 45%, of the study area exhibits a high susceptibility to landslides, as confirmed by the four machine learning models (Figure 6). Furthermore, the four susceptibility maps generated by machine learning show a uniform distribution of the different susceptibility levels, highlighting the study area's north-northeast regions as the most susceptible to landslides. Concurrently, the southern region is characterized as having a low susceptibility to landslides (Figure 5). The region exhibiting a notably elevated susceptibility to landslides is characterized by high precipitation, particularly in the north, along with contrasting relief. Flat and gentle slopes identify the low susceptible as gentle slopes. In summary, the landslide susceptibility maps produced in this study offer valuable resources for urban planners, government entities, and non-governmental organizations. These maps can aid in the effective management of human activities, resources, and infrastructure development within the region. The increased human activities in the region, combined with the landslide susceptibility maps generated, have the potential to enable local authorities to precisely identify areas at high risk and put in place early warning and preventive strategies for disaster mitigation.

Recent studies have noted a continuous increase in the frequency of extreme precipitation events in Mediterranean regions [START_REF] Fiori | Analysis and hindcast simulations of an extreme rainfall event in the Mediterranean area: The Genoa 2011 case[END_REF][START_REF] Tramblay | Future evolution of extreme precipitation in the Mediterranean[END_REF], translating into the more frequent occurrence of natural hazards, especially landslides and floods. The dominance of precipitation (DAAC between 16.72 and 14.66) over other factors in mapping landslideprone areas was expected since precipitation is regarded as the primary triggering factor for landslide movement and degradation, particularly in mountainous terrains (Argyriou et al., In nearly all machine learning models, after precipitation, the secondary factors influencing landslide susceptibility mapping include slope, distance to rivers, elevation, and distance faults (Figure 7). These four factors are associated with the geo-morphology characteristics of the study area, as the Rif region is known for its rugged topography, rough relief, and steep slopes. Geologically, it is a seismically active area with active regional faults and large shearing areas [START_REF] Poujol | Active tectonics of the Northern Rif (Morocco) from geomorphic and geochronological data[END_REF]. Such conditions made the drainage network flow mainly on fault corridors, deep valleys with steep slopes, following their orientation and adopting their structuring, leading to the formation of several types and patterns of drainage, especially trellis, rectangular, and barbed network patterns, which reflect an intense structural geological control on the drainage network [START_REF] Bouramtane | Drainage Network Patterns Determinism: A Comparison in Arid, Semi-Arid and Semi-Humid Area of Morocco Using Multifactorial Approach[END_REF]. As a result, the development of landslides frequently begins with the incision of rivers along the main faults and valleys, which defines the rhythm of the landscape's evolution. The slopes of the hills then respond to this incision by destabilizing the slopes' base and increasing the moisture content of the soil, leading to landslides, and risk is thus concentrated at the drainage networks (Figure 8) [START_REF] Clapuyt | Spatio-temporal dynamics of sediment transfer systems in landslide-prone Alpine catchments[END_REF]. Furthermore, this frequently observed spatial proximity of landslides to drainage networks is known as connectivity; it is crucial for large landslides in high-relief areas, and it is one of the distributional characteristics of landslides that determine their role as a natural hazard [START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF][START_REF] Roback | The size, distribution, and mobility of landslides caused by the 2015 Mw7.8 Gorkha earthquake, Nepal[END_REF].

The incorporation of the "distance to faults" factor among the prominent variables implies an association with the landslides identified in the study area, particularly those situated in close proximity to the major faults within the study area (Figure 8). This indicates the influence of tectonic activity on landslide risk, particularly the occurrence of co-seismic landslides. These co-seismic landslides are a significant secondary natural hazard in earthquake-prone mountainous regions, such as the Rif structural zone (Hovius et al., 2011; [START_REF] Roback | The size, distribution, and mobility of landslides caused by the 2015 Mw7.8 Gorkha earthquake, Nepal[END_REF]. Indeed, the Rif structural zone, particularly its eastern part covering the province of Taounate, experienced its last major earthquake of Mw 6.4 in 2004 and its most recent Mw 5.2 earthquake on 07/01/2023, with several Mw 4 to 4.9 earthquakes occurring throughout 2022. However, no studies have been conducted to date that directly link the study area's seismic activity to the occurrence of its landslides.

Performance evaluation of machine learning models for landslide susceptibility mapping

Several statistical methodologies have been recently employed to model landslide susceptibility in some limited areas of Taounate province in the Rif Mountains, including Ain aicha, Ourha, Oudka, and Taounate City. However, compared to conventional statistical analysis, machine learning methods demonstrate enhanced efficiency in addressing real-world challenges, such as landslide occurrences (B. Zhao et al., 2022). Machine learning methods have become increasingly valuable tools for decision-makers as they facilitate the identification of relationships between geological, topographical, climatic factors, human activities, and landslide events. Hence, comparative studies are needed to evaluate the performance of the models under the same conditions and make a fair judgment about their capabilities [START_REF] Rahmati | Land subsidence modelling using tree-based machine learning algorithms[END_REF]. In this present study, four machine learning models, namely SVM, CART, LR, and MDA, were compared and evaluated to map landslide susceptibility and identify the most significant causative factors. All four machine learning methods demonstrated excellent performance, highlighting the robust capabilities of machine learning modeling in the study area, with AUC scores ranging from 0.954 to 0.971 (Figure 4). The CART model demonstrated superior performance with an AUC value of 0.971. However, based on the evaluation metrics such as accuracy (89.92%, 83.01%), precision (81.74%), specificity (79.46%), and F1-measure (81.66), the SVM model performs at the highest level ( The CART model exhibited superior performance compared to all other models, as shown by the AUC and sensitivity indexes (AUC=0.971; sensitivity= 90.24%). This can be attributed to the capability of decision tree models to effectively address nonlinear problems and work on highdimensional data accurately (Razi & Athappilly, 2005). However, decision trees tend to overfit and do not generalize well to new data, and a slight change in the data can produce a very different decision tree [START_REF] Hong | Assessing landslide susceptibility based on hybrid Best-first decision tree with ensemble learning model[END_REF]. The results do not necessarily indicate that LR and MDA have the lowest predictive ability since they demonstrated high accuracy in AUC and some classification tasks. LR performed the best in terms of the F1-measure index and the second-highest AUC rate of 0.963. However, MDA performed best in terms of accuracy and specificity of training data (93.07%, 93.92%), as well as the sensitivity of test data (91.6%), which is considered a very good result for a linear model (Figure 4, Table 2). Moreover, the advantage of using the nonparametric MDA algorithm lies in its ability to perform effectively in the presence of complex associations [START_REF] Pourghasemi | Landslide susceptibility assessment and mapping using stateof-the art machine learning techniques[END_REF]. The LR and MDA models are sensitive to linear correlation of influencing factors, which are not adept at modeling complex nonlinear problems (C. [START_REF] Zhou | Landslide susceptibility modeling applying machine learning methods: A case study from Longju in the Three Gorges Reservoir area, China[END_REF]. However, the fact that the LR and MDA results are highly satisfactory and comparable to those of the SVM and CART models suggests the presence of a likely linear correlation between landslide occurrence and geo-environmental factors, which are likely the primary drivers of the landslide issue in the study region. Furthermore, using LR or MDA has a significant advantage in efficiently interpreting their findings. This is due to the standardized coefficients generated (Figure 9) by the LR model and the correlation factors (Figure 10) identified by the MDA model, which allows the identification of the processes that condition natural hazards through their interpretation of the LR model.

CONCLUSIONS

Identifying robust and efficient models to reduce errors in landslide susceptibility modeling and delimitation of landslide-prone areas is essential. The primary objective of this study is to evaluate four machine learning models (SVM, CART, LR, and MDA) for the purpose of mapping landslide-prone regions and gaining insights into the geo-environmental processes governing landslides in one of the most impacted areas in North Africa. 255 landslides were identified and collected with 15 landslide conditioning factors to develop the four machine-learning models. The dataset was divided into ten random divisions: 30% of the dataset was allocated for testing, while 70% was reserved for training. Six statistical measures, namely accuracy, precision, sensitivity (recall), specificity, F1-measure (or F1-score), and receiver operating characteristics (ROC), were used to compare the four models. The findings of the causal significance analysis indicate that landslides are primarily triggered by rainfall. Most of these occurrences are observed on slopes near rivers and faults, highlighting the significant influence of precipitation on landslides in the Rif mountains. The comparison results demonstrate that all four models have excellent and relatively similar performance (ROC curve > 0.954), but the SVM model performs best. The SVM model achieved the highest accuracy and F1-measure for the training data and the highest accuracy, precision, and specificity for the test data. Although the prediction accuracy of LR and MDA models is not as high as that of SVM, they outperformed the other models in some tasks and offer the advantage of providing more interpretable results regarding landslide hazard conditioning processes. Areas with a very high susceptibility to landslides make up nearly half of the study area, ranging from 40% to 45%. These areas are situated in the north and northwest, characterized by heavy precipitation and high relief. All four landslide models demonstrated strong performance in landslide susceptibility assessment, with the SVM model standing out as the top performer. Consequently, it is recommended for use in the creation of improved landslide susceptibility maps to enhance landslide hazard management. This study offers valuable insights for decision-makers in devising disaster risk reduction strategies in North Africa, particularly in regions sharing similar geo-environmental and climatic conditions. For future studies, it is essential to account for the distinct characteristics and influential factors specific to each type of landslide when crafting predictive models for their occurrence. Additionally, the combination of SVM with LR or MDA holds the potential to enhance model performance and offer a more comprehensible interpretation of landslide mechanisms.

Figure 1 .

 1 Figure 1. Location maps of the study area showing the location of landslides and non-landslides areas

Figure 2 .

 2 Figure 2. Methodological flowchart

Figure 3 .

 3 Figure 3. Landslide conditioning factors maps

  Lombardo et al., 2016; Ramos-Bernal et al., 2019). The jackknife test is employed to ascertain the individual contributions of each factor to the overall susceptibility model (Ramos-Bernal et al., 2019). The sensitivity of each of the fifteen factors in each predictive model has been assessed using the percentage decrease in overall accuracy (DAAC) (Bouramtane et al., 2022). 𝐹𝐹1 -𝑚𝑚𝑆𝑆𝐴𝐴𝑆𝑆𝐴𝐴𝐴𝐴𝑆𝑆 = 2 × 𝑇𝑇𝑇𝑇 + 𝐹𝐹𝑁𝑁 +

Figure 4 .

 4 Figure 4. The area under curve (AUC) related to the five employed models: support vector machine (SVM), classification and regression trees (CART), logistic regression (LR) and multivariate discriminant analysis (MDA)

Figure 5 .Figure 6 .

 56 Figure 5. Landslide susceptibility maps of Taounate province using the four models (a) CART, (b) SVM, (c) LR and (d) MDA

  2022; D'Ippolito et al., 2023; Tsunetaka, 2021). These results and findings are in line with the climatic characteristics of the study area, which includes both the Rifian region and the wider Mediterranean region. These areas are distinguished by substantial rainfall during the wet season (from November to February) and during summer, often in the form of storms (Senoussi et al., 1999).

Figure 7 .

 7 Figure 7. The contribution of the explanatory variables using DAAC for (a) CART, (b) SVM, (c) (LR) and (d) MDA models

Figure 8 .

 8 Figure 8. A map displaying landslides near faults and rivers: (a) near a fault and (b) near a stream

Figure 9 .Figure 10 .

 910 Figure 9. Standardized logistic regression coefficients of LR model

Table 1 . Details of landslide conditioning factors data

 1

	Cluster	Conditioning factors	Type	Data source
		Altitude	Continuous
		Slope	Continuous
		Aspect	Continuous
	Topographical	Curvature Profile curvature	Continuous Continuous	Digital Elevation SRTM -30 m
		Plan curvature	Continuous
		TPI	Continuous
		TWI	Continuous
	Geological	Lithology	Categorical	Geological data of Africa USGS + Geological map of Morocco 1/200,000
		Distance to faults	Continuous	Geological map of Morocco 1/200,000 + satellite imagery
	Hydrological	Precipitation Distance to Rivers	Continuous Continuous	CHIRPS -Satellite precipitation product (~ 5 km) Digital Elevation SRTM -30 m
		NDVI	Continuous	Sentinel-2 time series of 2021 -10 m
	Land covers	Land cover/Land use	Categorical	Esri Land cover/land use of 2021 -10 m
		Distance to roads	Continuous	Google Earth

Table 2 .

 2 The performance of four models in the terms of evaluation measures

		MDA		LR		SVM		CART
	Measures	Training	Test	Training	Test	Training	Test	Training	Test
	Precision	93.07%	70.12%	90.15%	81.62%	89.89%	81.74%	87.59%	71.43%
	Specificity	93.92%	61.77%	90.36%	78.06%	89.97%	79.46%	87.42%	65.82%
	Sensitivity	82.56%	91.60%	89.16%	90.29%	89.94%	89.24%	90.24%	77.69%
	Accuracy	88.24%	75.95%	89.69%	83.01%	89.92%	83.01%	89.10%	69.02%
	F1-measure	77.79%	60.84%	81.24%	79.12%	81.66%	72.48%	80%	54.67%

Table 2 )

 2 . SVM is, therefore, the most appropriate choice for mapping landslide susceptibility in mountainous regions. In this present study, our results align with earlier research conducted by Naceur et al., 2022, who employed models like Weight of Evidence (WoE), Radial Basis Function Network (RBFN), and SVM in the high Atlas of Zhou et al., 2018). It is suitable for predicting high-dimensional data, can also work on small data sets, and is robust when dealing with slight sample variations (Youssef & Pourghasemi, 2021). However, SVM is not an interpretable model, and it consumes much memory and time compared to other models (Marjanović et al., 2011).

	Morocco; Yousefi et al., 2020, conducted research
	in a mountainous region in Iran, utilizing General-
	ized Linear Model (GLM); Functional Discrimi-
	nant Analysis (FDA) and SVM; Pham et al., 2016,
	who applied Naïve Bayes (NB), Bayesian Network
	(BN), Naïve Bayes (NB), LR, SVM and Fisher's
	Linear Discriminant Analysis (FLDA) in India;
	Xie et al., 2021 conducted a study in China to as-
	sess landslide susceptibility using various methods,
	including SVM, and subsequently compared their
	performance. The SVM model achieved greater
	predictive accuracy due to its reliance on mini-
	mizing structural risk over traditional empirical
	risk minimization, aiming for a globally optimal
	solution for solving complex nonlinear problems
	(C.


Marilyn Keller

Marcell Krall

James Smith

Hans Clement

Alexander M Kerner

Andreas Gradischar

Ute Schäfer

Michael J Black

Annelie Weinberg

Sergi Pujades

Optimizing the 3D Plate Shape for Proximal Humerus Fractures

Keywords: Personalized plates, 3D printing, ex-vivo evaluation

To treat bone fractures, implant manufacturers produce 2D anatomically contoured plates. Unfortunately, existing plates only fit a limited segment of the population and/or require manual bending during surgery. Patient-specific implants would provide major benefits such as reducing surgery time and improving treatment outcomes but they are still rare in clinical practice. In this work, we propose a patient-specific design for the long helical 2D PHILOS (Proximal Humeral Internal Locking System) plate, used to treat humerus shaft fractures. Our method automatically creates a custom plate from a CT scan of a patient's bone. We start by designing an optimal plate on a template bone and, with an anatomy-aware registration method, we transfer this optimal design to any bone. In addition, for an arbitrary bone, our method assesses if a given plate is fit for surgery by automatically positioning it on the bone. We use this process to generate a compact set of plate shapes capable of fitting the bones within a given population. This plate set can be preprinted in advance and readily available, removing the fabrication time between the fracture occurrence and the surgery. Extensive experiments on ex-vivo arms and 3D-printed bones show that the generated plate shapes (personalized and plate-set) faithfully match the individual bone anatomy and are suitable for clinical practice.

Introduction

Bone shapes vary significantly across the world's population [START_REF] White | The human bone manual[END_REF], making the design of personalized medical plates for repairing fractures challenging. Yet patient-specific implants would provide major benefits for surgeons and patients: surgery time would be reduced and treatment outcomes improved. While some plates are malleable enough to be contoured during surgery to improve overall fitness, these are more prone to fatigue failure when compared to rigid precontoured plates. But anatomically pre-contoured plates only fit a subset of the population [START_REF] Goyal | Analysis of anatomic periarticular tibial plate fit on normal adults[END_REF][START_REF] Schmutz | Anatomical fitting of a plate shape directly derived from a 3D statistical bone model of the tibia[END_REF][START_REF] Kozic | Optimisation of orthopaedic implant design using statistical shape space analysis based on level sets[END_REF], causing malalignment at the fracture site [START_REF] Hwang | Mismatch of anatomically pre-shaped locking plate on asian femurs could lead to malalignment in the minimally invasive plating of distal femoral fractures: a cadaveric study[END_REF]. Even with the extra contouring step during the surgery, the plates do not fit all bone shapes [START_REF] Bou-Sleiman | Minimization of intraoperative shaping of orthopaedic fixation plates: a population-based design[END_REF].

Our work contributes several key steps to create personalized plates: (i) we first design an optimal PHILOS 3D plate shape for the humerus bone, fulfilling the clinical constraints given by a senior surgeon; (ii) we propose an anatomyaware transfer of the plate shape into any new bone, resulting in a personalized plate shape; (iii) we propose a method to position an arbitrary plate on an arbitrary bone according to clinical constraints, allowing surgical compatibility to be evaluated; (iv) we leverage the compatibility assessment to obtain a compact set of plates that accommodates a given population of bones; (v) we validate our methodology with extensive experiments on ex-vivo and 3D-printed bones, demonstrating the relevance of the designed personalized plate shapes for the clinical setting; (vi) we make available for research purposes the humerus and plate 3D models, as well as the plate extraction and fitting code † .

The creation of surgical plates consists of several tasks, which have been partially addressed by the existing state of the art. In Tab. 1 of Sup. Mat. we provide a schematic comparison of how our approach goes beyond existing methods.

Most methods start with a base plate or template that is designed for a specific bone surgery (tibia, humerus, clavicle, etc.) [START_REF] Goyal | Analysis of anatomic periarticular tibial plate fit on normal adults[END_REF][START_REF] Kozic | Optimisation of orthopaedic implant design using statistical shape space analysis based on level sets[END_REF][START_REF] Bou-Sleiman | Minimization of intraoperative shaping of orthopaedic fixation plates: a population-based design[END_REF][START_REF] Schmutz | Fit assessment of anatomic plates for the distal medial tibia[END_REF][START_REF] Schmutz | Fit optimisation of a distal medial tibia plate[END_REF][START_REF] Harith | Can we safely deform a plate to fit every bone? Population-based fit assessment and finite element deformation of a distal tibial plate[END_REF][START_REF] Carrillo | A time saver: Optimization approach for the fully automatic 3d planning of forearm osteotomies[END_REF][START_REF] Wu | Preliminary exploration of a quantitative assessment index for the matching performance of anatomical bone plates using computer[END_REF][START_REF] Tkany | New design process for anatomically enhanced osteosynthesis plates[END_REF][START_REF] Zenker | Quantifying osteosynthesis plate prominence-mathematical definitions and case study on a clavicle plate[END_REF]. One important step is how to position a given plate on a given bone. Existing approaches are either manual, semi-automatic [START_REF] Goyal | Analysis of anatomic periarticular tibial plate fit on normal adults[END_REF][START_REF] Schmutz | Anatomical fitting of a plate shape directly derived from a 3D statistical bone model of the tibia[END_REF][START_REF] Schmutz | Fit assessment of anatomic plates for the distal medial tibia[END_REF][START_REF] Wu | Preliminary exploration of a quantitative assessment index for the matching performance of anatomical bone plates using computer[END_REF], or fully automatic [START_REF] Kozic | Optimisation of orthopaedic implant design using statistical shape space analysis based on level sets[END_REF][START_REF] Harith | Can we safely deform a plate to fit every bone? Population-based fit assessment and finite element deformation of a distal tibial plate[END_REF][START_REF] Carrillo | A time saver: Optimization approach for the fully automatic 3d planning of forearm osteotomies[END_REF][START_REF] Tkany | New design process for anatomically enhanced osteosynthesis plates[END_REF][START_REF] Schulz | Evidence based development of a novel lateral fibula plate (variAx fibula) using a real CT bone data based optimization process during device development[END_REF]. Our fully automatic positioning strategy differs from the state of the art as it considers the surrounding anatomy; i.e. to minimize the risk of radial nerve damage, the plate is twisted around the bone so that the proximal end is fixed on the lateral side of the humeral head, and the distal end of the plate to the ventral surface of the humerus [START_REF] Gill | The spiral compression plate for proximal humeral shaft nonunion: a case report and description of a new technique[END_REF][START_REF] Da Silva | Comparing iatrogenic radial nerve lesions in humeral shaft fractures treated with helical or straight philos plates: a 10year retrospective cohort study of 62 cases[END_REF].

Once the plate is positioned, one needs to evaluate if the plate can be used in surgery. Some works compute plate-to-bone metrics [START_REF] Goyal | Analysis of anatomic periarticular tibial plate fit on normal adults[END_REF][START_REF] Bou-Sleiman | Minimization of intraoperative shaping of orthopaedic fixation plates: a population-based design[END_REF][START_REF] Zenker | Quantifying osteosynthesis plate prominence-mathematical definitions and case study on a clavicle plate[END_REF][START_REF] Petersik | A numeric approach for anatomic plate design[END_REF], but do not propose a binary decision criterion stating whether the plate is valid for surgery or not. We argue that this binary decision, also provided by existing works [3, 4, 7-9, 11, 12, 14], is important, as aggregate distance numbers can be misleading; e.g. one part of the plate could have a perfect fit but another part could make the plate not suitable for surgery. In our work, we define a three section-based criterion that takes into account the two plate-to-bone fixation regions and the middle plate section that transitions from one fixation region to the other. In addition, we perform a case study evaluation with an expert surgeon, in which the validity of the numerical criterion is confirmed.

Once the plate is positioned on a bone, and the fit evaluated, one can further consider the question of how to deform the plate to fulfill the fit criterion. For instance, [START_REF] Bou-Sleiman | Minimization of intraoperative shaping of orthopaedic fixation plates: a population-based design[END_REF] aims to develop a plate shape to reduce in-situ plate manipulation. They approximate the plate and target bone surfaces by planar sections and measure the bending necessary for each section of the plate to fit the bone surface. In contrast, our approach is similar to the one of [START_REF] Schmutz | Anatomical fitting of a plate shape directly derived from a 3D statistical bone model of the tibia[END_REF], in which the plate shape is not deformed after evaluating the fit, but rather directly extracted from the shape of an individual bone. We propose a plate transfer strategy that takes into account the matching of anatomic regions. One advantage over existing methods is that our approach simultaneously ensures both: (i) a personalized plate shape matching the bone's shape and (ii) its proper placement.

Prior work proposes a set of plate shapes that, together, can accommodate a range of patients [START_REF] Schmutz | Fit optimisation of a distal medial tibia plate[END_REF][START_REF] Harith | Can we safely deform a plate to fit every bone? Population-based fit assessment and finite element deformation of a distal tibial plate[END_REF]. This set is created by manually modifying an original plate to improve its plate-to-bone distance. We propose an automatic approach in which we first create many plates, and then, with a greedy algorithm, select the plate set that accommodates the most bones.

Method

Dataset. To create and evaluate our plate designs, we use a dataset consisting of 97 3D meshes of humerus bones, divided into two groups (A and B). Group A has 54 bones scanned with a FARO laser scanner (25 females and 29 males, 50% Black and 50% White, age range 17 and 45). Group B has 43 Computed Tomography (CT) scans of bones from autotomized body donors. The CT volumes were segmented and cleaned to reconstruct a mesh of the bone. Left-side humerus scans were mirrored along the z-axis to work with right-side humeri only.

Plate design. Given a bone, our goal is to automatically generate an optimal plate shape and determine its position on the bone. The optimality of our plate design is defined by an experienced surgeon, who established fixation points, areas to avoid, and plate-to-bone distance tolerances. These choices were validated by a second expert surgeon. A plate fulfilling these constraints on a bone is considered optimal for surgery. To obtain the plate design we 3D printed 7 bone meshes from our dataset with diverse shapes and asked the surgeon to annotate each bone with a marker. The surgeon indicated the bone areas that should be in contact with an ideal plate and the areas to be avoided (Fig. 1 left). From these annotations, we designed an ideal plate mesh P T contoured to a humerus template bone mesh T (Fig. 1 middle). This plate has similar dimensions as the actual INTEOS PROXIMAL HUMERAL PLATE 3.5. Following the surgeon's advise, the plate should not be in contact with the bone at the bone neck level, hence, an offset was added in this region. The fit criteria between a plate and a bone was defined such that the distance between them should not exceed 2mm at the fixation points, located at the head P h and the shaft P s , and should be less than 5mm in the other regions, with the exception of the plate neck, which must not be in direct contact with the bone (Fig. 1

right).

Bone registration for optimal plate transfer. Once the optimal plate is defined on the template bone, the goal is to register any new bone to the template, by preserving the relevant anatomic regions for the plate fixation. To do so, we adapt the registration technique designed for vertebrae [START_REF] Meng | Learning a Statistical Full Spine Model from Partial Observations[END_REF] to the humerus.

We first register each bone of set A by optimizing their Eq. 1, 2 and 3 from [START_REF] Meng | Learning a Statistical Full Spine Model from Partial Observations[END_REF]. Given these initial registrations we build a statistical shape model using principal component analysis (PCA). This model is parameterized with the shape vector β, and constructs a bone mesh T(β): an array of size N b × 3 containing the mesh vertex positions. Second, we register the whole dataset using the learned bone shape model by minimizing

E(β, t, r, F; S) = E p2m (S, r • (T(β) + F) + t) + λ β ||β|| 2 + λ L ∆(T(β)) • F (1)
where S denotes the scanned bone mesh, t and r are, respectively, the 3D rigid translation and rotation applied to the bone mesh vertices, and F is a 3D pervertex offset applied to each mesh vertex. E p2m (S, .) is the point-to-triangle distance between the scan vertices and the mesh triangles. The function ∆(T) is the mesh Laplacian operator, which is used to regularize the offsets F. We minimize Eq. ( 1) in three successive steps: with respect to (t, r), β, and F. For each scan S i we obtain its corresponding bone registration B ′ i = T(β) + F. Due to the elongated shape of the humeri and their axial rotation similarity, the registrations B ′ i can contain sliding, i.e. the same vertex does not correspond precisely to the same anatomic location on two different bones. To correct this, we perform a second pass of registration using smooth shells [START_REF] Eisenberger | Smooth shells: Multi-scale shape registration with functional maps[END_REF] to deform the initial template T to the surface of B ′ i and we obtain the final bone registrations B i matching the scan S i (see Fig. 2 left). Note that smooth shells [START_REF] Eisenberger | Smooth shells: Multi-scale shape registration with functional maps[END_REF] can not directly register T to the scans S i , and the intermediate meshes B ′ i are required. To validate the bone correspondences we define anatomic regions on the bone template T, and transfer them to all registrations B i on a per-vertex index basis. Fig. 2 shows that the annotated anatomic parts are faithfully preserved.

Plate extraction. As the registration process preserves anatomic regions, given a bone registration B i , we can extract a new optimal plate P i . For each vertex of the designed plate P T we compute its offset to the template bone T. By applying these offsets to the new bone B i , we obtain the personalized plate P i .

Positioning a plate on a bone. Now that we can create personalized plates, we want to study if a created plate can be used for surgery on another bone. For that, given an arbitrary plate P i and bone registration B j , we need to position the plate and determine if it fits. We automate the plate positioning by minimizing the distance between the plate and bone fixation points. We start by computing an initial 3D rigid transformation r 0 , t 0 by optimizing Eq. 2 in Sup. Mat. and then obtaining the final positioning r f , t f by optimizing Eq. 3 in Sup. Mat. These equations enforce the fixation areas of the plate P i h to match B j h and P i s to match B j s (see Fig. 1 right), while avoiding plate-to-bone inter-penetrations. Once the plate is positioned on the bone, we evaluate the binary fit criteria (see Fig. 1 right) to conclude whether the plate shape is fit for surgery or not.

Plate shapes set.

A single plate design can not accommodate the whole population due to its morphological variance [START_REF] White | The human bone manual[END_REF]. Thus we propose to build a set of plate shapes with a greedy algorithm. The algorithm is summarized in pseudocode in Alg. 1 of Sup. Mat. Given a bone dataset, we start by creating as many personalized plates as bones. Then, for each plate we use the previous optimization to determine how many bones it accommodates, i.e how well a plate shape generalizes to different bones. We then select the plate that fits most bones and remove the plate and the fitted bones from the current sets. We iterate until no bones are left. The result is an ordered plate-set that accommodates the input bone dataset.

Experiments

Numerical evaluations

Bone registration accuracy. We validate that our registrations B i accurately match the scans S i by computing the mean distance (MD) between each registration to the closest vertex in the scan. We obtain a MD of 0.08 mm (std=0.04) for group A, and a MD of 0.27 mm (std=0.55) for group B. For both sets, the registered bones match the scans with sub-millimeter accuracy. In addition, all distances are less than 1mm on the surface where the plate is positioned (Fig. 3 left). Plate set coverage. One application of our design is that it can generate a plate set that could be preprinted in the hospital and be readily available for immediate use. Using our greedy algorithm on our bone dataset, a set P N =5 S made of the first 5 plates accommodates already 51.04% of the bones and P N =10 S accommodates 73.96% (Fig. 3 middle). In Sec. 4 we discuss how this coverage could be improved.

Ex-vivo evaluations

Comparison with state-of-the-art plates. We CT-scanned an isolated bone (not included in the original dataset) and asked a surgeon to manually contour a state-of-the-art plate, as done in clinical practice. We compare it to the custom plate P C and the best set plate P S by computing the plate-to-bone distances (Fig. 3 right). The proposed plates are clearly closer to the bone than the handcontoured plate. This proximity is known to be beneficial for bone recovery [START_REF] Ahmad | Biomechanical testing of the locking compression plate: when does the distance between[END_REF]. While P S is less accurate than the custom plate -and shorter, as it was generated from a shorter bone -the surgeon considered the plate and its placement as suitable. This experiment reveals the difficulty of closely fitting the bone by manually bending a plate and the benefit of the proposed design.

Ex-vivo surgery. We performed an ex-vivo experiment on a cadaveric arm, mimicking an actual minimally-invasive surgical operation setting. One goal was to test whether the designed 3D plates can be properly inserted along the bone under the muscles. We CT-scanned 3 cadaveric arms (not included in the original dataset), reconstructed the bone scans and registered them to obtain B 1 , B 2 and B 3 . For each bone we generated and 3D printed their custom plate P 1,2,3 C and the best fitting set plate P 1,2,3 S . The plates were drilled and coated with metallic paint to be visible on the CT scan. During the experiment, the plates Fig. 4: Left: Color-coded plate-to-bone distances (blue: 0 mm; red: ≥ 5 mm) for the set plate and custom plate on bones B 1 , B 2 , and B 3 . Right: Placement of the 3D printed plates P S on the corresponding 3D printed bone: green: algorithm's placement; orange: the surgeon's placement.

were inserted in the cadaver arm, as it would be done in a standard operation, and fixed to the humerus with screws. For the 3 bones, the surgeon estimated that the fit of both the custom and set plates was satisfactory. He noted that the surface roughness of the 3D print and metallic coating was too high, making the plate hard to slide across the humerus into place.

We CT-scanned each plate screwed to the bone and used the scans to asses the fit quality. Fig. 4 left shows the plate-to-bone distances for each bone. The plate P 1 C fits closely on the proximal side near the sulcus bicipitalis. For P 2 C and P 3

C the fits are not as good as P 1 C . This was credited to a handling/alignment issue raised by the surgeon. Consequently, greater deviation in alignment is seen on the distal end of the bone on the lateral side. There is also a greater distance between the ventral part of the plate and the sulcus bicipitalis. The plate is too lateral on the humerus. In each case, the surgeon found the custom plate better. The limit in terms of fit was mainly attributed to the operative process rather than to the plate shape itself. The surgeon was satisfied with the 6 fits and noted that the 3D-printed plates were less stiff than the commercial plates he usually uses, making them harder to manipulate.

Evaluations with 3D-printed bones

To further evaluate our design, we generated 7 bones {B i } i∈ [4;10] . B 7 is the mean bone shape and the 6 others are generated by uniformly sampling the bone PCA space B (β ∈ {-2, 2}) in the first 3 PCA dimensions. For each bone, we 3D printed its custom plate, P C , and the best set plate, P S , from the plate-set.

First, we presented pairs of bones and plates to the surgeon and asked whether they were fit for surgical use. Then, we asked them to position the plates on the 3D bones to evaluate our positioning strategy. Last, to evaluate the adequacy of a plate set, we presented one bone to the surgeon and asked them to choose the best plate from the set.

Ready-for-surgery evaluation. For each bone B i , the surgeon considered both plates as candidates and evaluated if their shape was suitable for surgery. The surgeon stated that 100% of the plates were fit for surgery and that the custom plates P C were always clearly better fitting than the plates P S from the set.

Positioning evaluation. To evaluate the plate positioning computed by our algorithm, we compare it to the surgeon's positioning. We asked the surgeon to position the set plates P S on each bone, secured them with blue tack and scanned the whole with a FARO hand scanner. Fig. 4 right shows the overlay of both positions of the used plate for each bone. For most bones, the placement of the plate by the surgeon and the algorithm were very similar. The differences can be explained as there are multiple plate positions that fit. Both placements (orange and green), even in the most dissimilar case (first bone in Fig. 4 right) were considered correct by the surgeon.

Choosing from the plate set. Our algorithm iteratively positions the plates P i S from i = 0 to i = N and picks the best-fitting plate, so we asked the surgeon to perform a similar trial-and-error task to evaluate the relevance of our plate set.

Out of the 7 bones, the surgeon only chose the same plate as our algorithm for two of them (B 7 and B 10 ). For the other cases, the surgeon compared his selection to the algorithms' choice and confirmed that both plates were fit for surgery: different plates with slightly different placements can work for the same bone. Most interestingly, some plates selected by the surgeon did not fit the bone according to the strict numeric fit criteria. We discuss this finding in the next section. The surgeon noted that the plates of our set were always preferred over the state-of-the-art plates.

Conclusions, take-aways and future work

Conclusion. We propose to automatically generate a custom humerus plate that specifically matches the 3D shape of a bone and specifically meets the requirements of a surgeon for minimally invasive surgery. We also generate a set of plate shapes that accommodates a given bone population. We extensively evaluate our approach on cadaverous arms and 3D printed bones.

Experimental takeaways. All the proposed plates (individual and set plate) match the shape of the bones, while the individual ones are considered the best by the surgeon. Furthermore, our experiments show the importance of evaluating on 3D-printed bones and ex-vivo arms. With the former, the surgeon can easily assess the plate-to-bone fit. With the latter, the insertion procedure can affect the plate placement. Moreover, our results on 3D printed bones argue for a relaxation of the theoretical fitting constraints: the plate-to-bone distance could be higher in some areas while still obtaining a proper surgical fit. Loosening the tolerances would allow more bones to be fit with the same plates and potentially the number of plates required in the plate set could be further reduced.

Future work. To reach actual clinical use, two elements need consideration: (i) the current method takes as input the shape of a healthy bone, whereas patients have a fracture. One could perform a CT scan of the other arm and generate the symmetric bone or leverage methods that reconstruct a full bone from partial observations [21]. (ii) Our work focuses on the geometric design of the plate but does not consider the physical properties of the plate, such as thickness and stiffness. Future work should optimize these physical properties for the stresses that the plate must endure, e.g. using finite element methods.

Fig. 1 :

 1 Fig. 1: Left: surgeon annotations showing where the plate should be positioned and areas to be avoided. Middle: the designed plate on the template bone. Right: Plate-on-bone fit criteria. Colors show the vertices defined as the fixation area on the plate (Blue: P h ; Red: P s ) and on the bone (Pale blue: B h ; Pale red: B s ).

Fig. 2 :

 2 Fig. 2: Left: Registered bone B i (green) and superposition to the input scan S i (pink). Right: Registrations B i to multiple scans. Despite the different individual bends and twists of each bone, the colored anatomic regions are consistent.

Fig. 3 :

 3 Fig. 3: Left: Per vertex max distance between the registrations B i and scans S for bone set A (left) and B (right). Blue (0 mm) to red (1 mm). Max distances are higher for bone set B because of their lower scan quality (segmentation artifacts). Middle: Accommodation percentage as a function of the number of plates. Right: Plate-to-bone distance for three plates: surgeon's hand bent, best from the plate set, and custom plate. Blue (0mm) to red (5mm or more).

  bone and implant significantly reduce construct stability? Injury 38(3), 358-364 (2007) 21. Meng, D., Keller, M., Boyer, E., Black, M., Pujades, S.: Learning a statistical full spine model from partial observations. In: International Workshop on Shape in Medical Imaging. pp. 122-133. Springer (2020)

† https://humerusplate.is.tue.mpg.de

Acknowledgements

This work was supported by the project CAMed (COMET K-Project 871132) which is funded by the Austrian Federal Ministry of Transport, Innovation and Technology (BMVIT) and the Austrian Federal Ministry for Digital and Economic Affairs (BMDW) and the Styrian Business Promotion Agency (SFG). Michael J. Black (MJB) has received research gift funds from Adobe, Intel, Nvidia, Meta/Facebook, and Amazon. MJB has financial interests in Amazon, Datagen Technologies, and Meshcapade GmbH. MJB's research was performed solely at MPI.

For the scanned humerus dataset, we acknowledge Dr. Brian Corner, Dr. Jeff Hudson who did the scanning, and the Cleveland Museum of Natural History Physical Anthropology Dept for access to the Hamann-Todd human skeletal collection. We also thank Anja Gieringer for the bone CT scans and Karoline Seibert at Hofer GmbH & Co KG, Fürstenfeld, Austria for the 3D printing of bones and plates.


A NEURAL-NETWORK-BASED SYNTHETIC BOUNDARY CONDITION FOR THE SIMULATION OF CHANNEL TURBULENCE

Heesik Yoo 1 , Bérengère Podvin * 1 , Yann Fraigneau 1 , and Alexandre Allauzen 1

1 Université Paris-Saclay, CNRS, LIMSI Summary We describe a method combining Proper Orthogonal Decomposition (POD) and an artificial recurrent neural network to implement an inlet boundary condition for a turbulent channel flow. The boundary condition is reconstructed from POD modes, the amplitudes of which are predicted using a neural network (LSTM). We first assess the quality of the LSTM prediction in POD space. The reconstructed velocity is then used as an inlet boundary condition in a spatially developping channel. Statistics are compared with the reference and are found to be in relatively good agreement.

MOTIVATION

The generation of turbulent inlet boundary conditions is an important question for the simulation of flows around aerodynamic immersed bodies or a spatially developing boundary layer. Fully developed turbulent channel flow can be obtained using periodic boundary conditions, but in general turbulent flows such as boundary layers develop in the spatial direction and more effort is required to construct an inlet velocity condition.

Pamiès et al. [START_REF] Pamies | Generation of synthetic turbulent inflow data for large-eddy simulation of spatially evolving wall-bounded flows[END_REF] used a combination of structures to construct the inlet condition -synthetic eddy method. Proper Orthogonal Decomposition (POD) provides a natural reconstruction framework as the velocity field can be approximated as a linear combination of POD spatial modes, the amplitudes of which vary in time. The question of determining inlet boundary conditions thus boils down to predicting the temporal evolution of a relatively limited number of amplitudes. In this framework the issue is to predict the behavior of the amplitudes of the POD spatial modes. Druault et al. [START_REF] Druault | Generation of three-dimensional turbulent inlet conditions for large-eddy simulation[END_REF] combined linear stochastic estimation with POD to construct an inlet boundary condition for a spatially developing mixing layer. Perret et al. [START_REF] Perret | Turbulent inlet boundary conditions for large-eddy simulation based on low-order empirical model[END_REF] used experimental data to construct a low-order model for the amplitudes of the POD modes and build an inlet condition for large-eddy-simulation (LES) of a mixing layer. In the present work, we predict POD amplitudes using a special type of neural network called a long-term short-term memory network (LSTM) [START_REF] Bucci | Control-oriented model learning with a recurrent neural network[END_REF]. We first describe the method and show how the network is able to predict the POD amplitudes of the velocity field in a turbulent channel flow. The method is then implemented to reconstruct the velocity field on a plane which is used as an inlet boundary condition for the numerical simulation of a turbulent channel flow. Turbulent statistics in the channel are then compared with those of a reference one.

RECONSTRUCTION METHOD Description

We consider the numerical simulation of an incompressible turbulent channel flow characterized by a bulk velocity U and a half-height h. The Reynolds number based on channel half-height h and friction velocity u τ is Re τ = 194, which corresponds to the reference channel. Application of Proper Orthogonal Decomposition to the velocity field yields

u(x, t) = n a n (t)φ n (x) (1)
The spatial eigenfunctions φ n are determined a priori from a reference simulation. We consider 300 snapshots separated by a time scale of δt = h/U , which can be reorganized into a hierarchy of 300 modes.

Constructing an inlet condition at all times consists in predicting the amplitude of these modes. In order to this we use a LSTM. The network is trained on a subset of the snapshots and validated on the remaining subset. The goal of the LSTM is to predict the next values of the amplitudes at t + δt, given the values at t. It is a recursive procedure, so that the amplitudes can be predicted over a long period of time.

Once the predicted values obtained at regular intervals of δt, the amplitudes can be interpolated in the simulation at any instant within the range spanned by the prediction.

Results

In the results shown below, the LSTM was trained on 180 snapshots. The remaining 120 snapshots constitute the validation basis. The LSTM was iterated over 640-180=480 snapshots -120 of which can be compared with the reference and 360 of which correspond to entirely new times. We consider only the 100 most energetic POD modes, which capture about 85% of the energy.

We first compare the amplitude of the POD modes predicted with the LSTM to those the reference simulation. Figure 1 represents respectively the energy content and frequency spectrum of the amplitudes (frequencies are expressed in wall units, denoted with a +). Overall the LSTM predicts relatively well the characteristics of the modes.

INLET BOUNDARY CONDITION

In this section we simulate a spatially developing turbulent channel using as an inlet condition the velocity field reconstructed using the POD eigenfunctions and the amplitudes predicted by the LSTM. A convective boundary condition is used at the outlet of the channel. The dimensions of the channel are L x = 2π, L -z = 4π, L y = 1 in respectively the streamwise (x), wall-normal (y) and spanwise (z) direction. The statistics are computed over 24 time units based on channel half-height and bulk velocity in cross-planes at the following streamwise locations: the inlet, the mid-channel x = π and the outlet. They are compared with those of a reference simulation at the same Reynolds number with periodic boundary conditions in the streamwise direction.

Figure 2 shows the turbulent intensities (or normal components of the Reynolds stresses). The normal stresses components of the reconstructed velocity condition are slighty under-predicted, which corresponds to the fact that only 85% of the kinetic energy is captured in the reconstruction. However at the outlet of the channel, the stresses are in much better agreement with the reference.

CONCLUSIONS

A new method combining Proper Orthogonal Decomposition and neural approaches was developed in order to build an inlet condition for an turbulent channel flow. The boundary condition was reconstructed using POD modes determined a priori. The amplitudes of the POD modes were predicted using a LSTM network. The characteristics of the amplitudes were satisfactorily captured by the LSTM. The statistics of the channel with the synthetic inlet boundary condition were found to match relatively well those of a turbulent channel flow with standard periodic boundary conditions.

Figure 1 :

 1 Figure 1: Left: Energy of the POD amplitudes as a function of POD mode number for the reference data and the prediction; Middle: Frequency spectrum of the reference POD amplitudes as a function of POD mode number; Right: Frequency spectrum of the predicted POD amplitudes.

Figure 2 :

 2 Figure 2: Turbulent intensities averaged over time and the spanwise direction at different streamwise locations as a function of the wall distance y (expressed in wall units +); top: channel inlet; bottom: channel outlet.