Nils Gentschen Felde
email: felde@nm.ifi.lmu.de

Mabel Fraume

Matti Heikkurinen

Dieter Kranzlmüller
email: dieter.kranzlmueller@lrz.de

Julio Serje
email: serje@un.org

Nils Gentschen Felde

Mabel Cristina

Marulanda Fraume
email: marulandafraume@un.org

UNISDR Global Assessment Report -Current and Emerging Data and Compute Challenges

Keywords: Hazard and loss modelling, probabilistic modelling, disaster risk, distributed data management

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

Introduction

This paper describes the data and compute challenges related to the production of the biennial Global Assessment Reports (GAR) [START_REF]UNISDR: Global Assessment Report on Disaster Risk Reduction (GAR)[END_REF], key documents providing highlevel overviews of the status of the disaster risk reduction activities on the global level. The production of the data these documents are based on -hazard scenarios, exposure information and vulnerability modelling -is performed by distributed collaborations and coordinated by UN Office for Disaster Risk Reduction (UNISDR). The risk calculation that provides estimates of the likely annual losses due to natural disasters represents a time-critical challenge, both in terms of organising the necessary computational processes to produce and verify the results, and in terms of managing the data sets in a consistent way. In parallel to the GAR-specific analysis, making the modelling data available as an open data service could support numerous additional research activities. These issues and goals form the context of the ongoing collaboration between UNISDR, LMU and LRZ.

While the type of large-scale disaster risk modelling GAR represents is most likely unique, we can see similarities with overall process and organisation as well as technical approaches with certain initiatives in other research domains. For example, the overall organisational structure resembles the approach used by global High-Energy Physics (HEP) collaborations. However, managing the complex, interdependent evolution of the interfaces (both physical and software ones) between thousands of components in a typical HEP project has necessitating developing relatively rigid and formalised organisational processes. This makes most of the tools developed to support HEP collaborations (ranging from document management systems -such as EDMS [2] -to global data/compute systems such as Worldwide LHC Computing Grid [START_REF]Worldwide LHC Computing Grid[END_REF]) developed for HEP collaborations not optimal for the GAR process.

Perhaps the closest analogue can be found from the earth observation domain. The Group on Earth Observations (GEO) has launched the GEO-DARMA [START_REF]GEO-DARMA = Data Access for Risk Management, Group on Earth Observations[END_REF] initiative with a goal of bringing earth observation data into disaster risk management. The effort builds on earlier initiatives focusing on specific hazards (floods, volcanoes) and aims to extend the focus from supporting the immediate, acute response to supporting preparedness and risk reduction. However, at the time of writing the initiative is still in its early stages.

2

Global Assessment Report -GAR

Background

UNISDR was established in 1999 and its role in the UN system is to serve as the focal point for disaster risk reduction activities to ensure coordination and synergies between UN organisations, and regional and national activities. The two major UN policy documents bringing all this guidance together into top-level policy documents are the 2005 Hyogo Framework for Action [START_REF]UN-ISDR: Hyogo Framework for Action 2005-2015: Building Resilience of Nations and Communities to Disasters. United Nations -International Strategy for Disaster Reduction[END_REF] and the 2015 Sendai Framework for Disaster Risk Reduction [START_REF]United Nations General Assembly: Sendai Framework for Disaster Risk Reduction 2015 -2030[END_REF]. The GAR process played a key role in implementing the Hyogo Framework and in the preparations of the Sendai Framework. It provided concrete data and examples of how the policies implemented (or to be implemented), changing natural conditions (e.g. climate change), population movements and major infrastructure projects are influencing the likely consequences (lives lost and direct economic losses) of natural hazards.

On the abstract level, the Global Risk Assessment processes of GAR 13 and GAR 15 are based on building three datasets used to generate risk metrics such as the "Loss Exceedance Curve" (LEC), the "Average Annual Loss" (AAL) and the Probable Maximum Loss (PML) plot the for the different hazards considered for each country. The data necessary for performing the global risk assessment of the GAR15report are:

• Hazard data, consisting of groups of simulated scenarios for each of the natural hazards for each of the natural hazard (earthquake, tsunami, riverine flood, cyclonic wind, storm surge and so on) used in the analysis. Each set of simulated scenarios must comply with the certain key requirements, such as being mutually exclusive, collectively exhaustive and having an annual frequency of occurrence associated with them. • Exposure data, describing each exposed asset with a set of attributes such as their geographical location, structural characteristics, construction material type, economic value (among others). • Vulnerability data, characterising the exposed asset with a set of attributes describing their relevant characteristics that determine how sensitive they are to different hazards at different intensities. Current GAR information linked to other datasets allow further evaluations and analysis. For example, risk associated to hydrometeorological hazards are strongly influenced by climate change, hence the IPCC data IPCC [START_REF]Intergovernmental Panel on Climate Change[END_REF] data and reports used as an input for simulation of new hazard scenarios considering climate change. Similarly, the exposure data used includes contributions collected using crowdsourcing approach based on OpenStreetmap [START_REF]OpenStreetmap Project[END_REF] and the vulnerability data can be improved by counting with better information of exposed assets. Additional data sources, such as OECD macroeconomic data, are used to conceptualise the disaster risk metrics. These dependencies are a partial rationale for the relatively frequent releases of the GAR.

The key output of the analysis are the country-level summaries presented in the Global Assessment report (a typical view of the online version is presented in Figure 1). This distils a complex and multi-faceted analysis into a summary with few key indicators that are suitable for steering policy-level decision making in the UN member states.

As a result, the environmental modelling behind GAR can have a major societal impact. While not legally binding, the GAR recommendations have an impact on national legislation and e.g. zoning decisions -both areas that unavoidably have an impact on economic situations and prospects of both public and private sector entities. In the long-term it has been shown that most of the investments in risk reduction are "profitable" in the sense that investments will eventually prevent direct economic losses that would have been several times higher than the money spent on protection (as an example, analysis of government-funded flood defence schemes in the UK showed an average benefit-cost ratio or 9.5:1 [START_REF] Fraume | Strategic Flood Risk Management[END_REF]). However, statistically major natural hazards usually have several decades between each occurrence, complicating the benefit-cost calculations. Justifying the immediate costs (loss of revenue or increased tax burden) by increased resilience in situations that statistically occur e.g. once every 50 years can be politically challenging.

For reasons outlined above, the results of the modelling can be expected to be under more scrutiny than a typical peer-review process for an academic publication. This pressure is further increased by the fact that the scale of the problem necessitates limiting the granularity of the analysis from what would be possible on local or regional analysis. Thus, the modelling software, loss calculation process and all the related data management practices need to be monitored very carefully.

GAR contents

As outlined earlier, the GAR analysis process needs to take into an account changes in the hazards themselves (e.g. increased frequency of extreme weather events), the developments in the distribution of population and infrastructure (e.g. urbanisation) and the impact of the policies implemented so far (e.g. changes in building codes, flood barriers, early warning systems and so on). Based on input data from UN member states and from other sources (such as WMO and IPCC), a large-scale simulation effort will produce tens of thousands to millions of hazard-specific scenarios (e.g. flood, seismic, tropical storms etc.) describing the location and intensity of the natural hazard. The hazard-specific scenarios are put together with the country-specific exposure (geographical location and physical attributes) and their associated vulnerability corresponding to each hazard in order to perform probabilistic calculation of risk (likely losses).

The resulting risk metrics such as the AAL and PML for different return periods are put into context, e.g. by comparing the average direct economic losses to key macroeconomic indicators (see Figures 1 and2). This can be used to gain a quick overview of the risks on the global scale (see Figure 3). An interactive viewer [START_REF]Global Assessment Report on Disaster Risk Reduction 2015, Risk Data Platform CAPRAViewer[END_REF] is also available online. The implementation of the Sendai framework for disaster risk reduction will increase the interest on both the GAR reports themselves, as well as the underlying data. The new framework for disaster risk reduction calls for more ambitious monitoring and modelling of risk. The plans to move into a "on demand" approach for risk modelling -as well as related policy developments that e.g. call for taking disaster risk into account when planning any investments -increase the demands on the assessment process. In addition, they make the provision of direct access to hazard, exposure and vulnerability data an important tool for increasing synergies between these activities. Fig. 3. A global summary of the proportional impact of the AAL to national economies

3

The GAR compute challenge

The GAR compute challenge can be broken into two parts: generating the hazard scenarios and combining them with the exposure database and the associated vulnerability functions to calculate the potential losses. The main dependencies between these steps are presented in Figure 4. The tools used in this process are heterogeneous, usually developed independently by the teams who are responsible for specific subtasks. Thus, the mode of operations is both globally distributed and very heterogeneous. On the technical level the CAPRA-GIS [11] toolkit plays an important integrating role: it provides the common formats for presenting hazard scenarios, as well as providing foundations for the overall loss calculation process. Hence the computational challenges have a clear interface between them.

Fig. 4. Organisation of the teams involved in GAR analysis

Production of hazard scenario files

The hazard scenarios are produced using a wide variety of methods, ranging from models running on powerful desktop computers to ones using computer clusters to perform the work. The basic probabilistic process is similar: the hazards scenarios are produced independent from each other and can thus be considered trivially parallelisable processes. The implementation tools range from software developed completely inhouse to models that run on platforms such as Matlab. However, despite the considerable computing resources needed, this step is rarely a time critical issue, as the schedule for producing the data is known well in advance and due to lack of dependencies between individual calculations any additional resources that can be brought in will speed up the process.

This situation may well change if the anticipated move to a more dynamic GAR process will be extended to the production of scenario files. For example, the overview document describing the latest approach to the production of the flood model data [START_REF] Rudari | Improvement of the Global Flood Model for the GAR 2015[END_REF] lists 19 different external datasets used to initiate and fine-tune models, which represents a part of the process that is inherently serial in nature, representing an execution step that will not benefit from additional computing resources.

Loss calculation

The performance of the loss calculation step tends to be the most time critical part of the Global Risk Assessment process -already with the 2-year publication cycle. Any unanticipated delay in the generation of the hazard scenarios will delay the overall loss calculation. Furthermore, due to the importance of presenting the data on per-country basis, summarising the loss estimates of large countries (such as China) can take a very long time.

Development of the CAPRA framework needs also to ensure that any new versions of the software produce same results as the original reference software versions. Hence, short-term developments tend to be incremental in nature. However, this incremental approach may face challenges with increased resolution (a uniform 1 km x 1 km grid instead of the current approach using resolutions ranging from 1 km x 1 km to 30 km x 30 km, depending on the hazard and geography of the region) and especially more dynamic production schedule of the GAR.

At the moment, the risk calculations are performed using two different versions of the CAPRA software: the original, single-threaded version (implemented in Visual Basic and available as part of the overall CAPRA-GIS package [11]) is used for most of the risks, as it is sufficiently powerful for risks where high-resolution modelling is not needed (e.g. modelling of droughts). A version of the risk calculation engine that has been parallelised and ported to Java (by the UNISDR team in collaboration with LMU and LRZ) is currently the reference implementation for flood modelling and is planned to be taken into use for the assessment of earthquakes in the near future. The Java version is already capable of exploiting shared memory parallelism and achieving close to a factor 50 speed up compared to single-threaded version. It is likely that further, incremental development of this version of the software will be sufficient to cope with the increase in resolution in the major publications (with the 2-year production cycle). However, additional optimisation is likely needed for the on-demand production of the reports (some potential approaches are discussed in chapter 6 of this paper).

4

The GAR data challenge

The current GAR data challenge can be broken down into two main phases: the hazard scenario stage and the data transfer to UNISDR. The vulnerability and exposure datasets are very small in comparison, and produced in more centralised manner. In either of the cases the amount of data is not a major problem per se, but rather the latencies introduced by the data production and transfer. While the hazard scenario development may in some cases need considerable resources, even the largest global hazard dataset will be of the order of few Terabytes. The situation might change slightly in the future due to higher resolutions used in the risk modelling. However, as the information is transmitted in compression format move from "5 km x 5 km" to "1 km x 1 km" resolution will not mean that the amount of data would automatically grow by factor of 25.

A bigger issue is moving the data to the loss calculation team in an efficient and coherent manner. While the current datasets are not excessively large (of the order of few Terabytes), limitations of the network infrastructure available to some of the partners in the collaboration necessitates moving the data by sending physical hard drives through postal or courier services. This approach is unlikely to cause insurmountable problems in the near future, as long as the resolution of the analysis will not be increased dramatically. It is likely that the capacities of commodity hard drives will grow at sufficient rate to match the increase of the data volume.

However, the approach is not without its issues: copying the data from the original storage system to a transient media may introduce issues with consistency, especially if there will be additional versions of the data that complement the biannual GAR process. Developing processes and metadata approaches to handle these issues are relatively straightforward to manage in setting where the data processing is done by a relatively small, established collaboration. However, turning the hazard, vulnerability and exposure datasets into open data products and services used by a broader research community will likely bring up additional issues that GAR collaboration needs to address e.g. through additional documentation or training activities.

Impact of the on-demand process

Traditionally the GAR process has required storing only a complete, global dataset for each of the biannual publications. This means that even considerable increase in resolution will not increase storage capacity needed beyond what is possible to handle using commodity solutions. The two-year cycle will also create a framework for managing versions of the datasets in a very intuitive manner: even a directory structure that is based on the production year of GAR is sufficient in most cases.

The on-demand process will create additional challenges. When parts of the analysis will complement the main GAR dataset with updated information related to a country or a region, the consistency of the data management and the ability of the metadata system to maintain link between a publication and the corresponding dataset will need to be reassessed. With the increased resolution it is also possible that -if these additional model executions and subsequent versioning of the result datasets are a frequent occurrence -the size of the overall dataset grows to a level that makes moving it challenging. This may create a need to re-examine the current distributed computing model, as in more and more cases performing the computation in the same computing centre that holds the relevant datasets are stored will be advantageous or even necessary from the performance point of view. Open data prototype and pilot

The motivation for investigating the feasibility of an open data approach are manifold, ranging from principles related to transparency of the GAR process to catalysing and supporting open innovation ecosystem that could also uncover novel approaches to disaster risk reduction. The data is already shared on request and used actively by third parties (e.g. by insurance companies to support their internal risk assessment processes). However, providing potential new users instant access is seen as a key method for removing barriers to new research and innovation activities.

The technical approach chosen for the open data pilot can be characterised as a "Minimal Viable Product" approach. The starting point is simple: a download portal addressing the key requirements of the GAR team and the external parties participating in the "beta phase" of the open data pilot:

• Download functionality: browse and download individual files or directories • Mechanisms for branding and ensuring that users are aware of the licensing issues and key disclaimers related to interpretation of the data. • Basic mechanisms for linking metadata to the data itself • To be considered: upload/updates through the web interface -maintain consistency between different copies and with metadata

The first versions of the portal were developed at LMU in fall 2016, and will be used to refine the requirements of the production version that may eventually become part of the formal LRZ service portfolio. Once sufficiently mature, the solution will most likely be merged with the current, download functionality [START_REF]UNISDR, Risk Data and Software Download Facility[END_REF] that is included in online version of the GAR, providing access to a limited subset of the information. The immediate needs for further development are mainly related to supporting archive file formats that are commonly used by the intended user community (possibly zip or rar formats in addition to tar.gz). In the longer term, aligning the data download functionality with the online visualisation tools and the Risk Atlas will also need to be considered.

Outside the purely technical issues, the open data pilot may bring up new requirements in terms of engaging with the community. Making the data available in a way that decouples it from the interpretations made by the experts (either in GAR or by the groups who developed the hazard models) may bring up situations where the interpretation -or misinterpretation -of the data outside the strict UNISDR scope should at least be brought to the attention of the original experts. Thus, supporting functionality such as social media integration may be deemed a useful complement to the download functionality itself.

Future directions

The open data portal will most likely play a role in the launch of the 2017 version of the Risk Atlas that provides a consistent geographical view of the disaster risk information, allowing business, investors and international organisations assess and compare the risks and resilience to the occurrence of natural hazard events in different countries more intuitively. Linking the GAR data organically to such an overview document will most likely increase interest in both. Thus, determining methods to discover and support new use cases emerging from the use of Risk Atlas and GAR data will proceed in parallel with the implementation of the first version of the open data service.

We foresee that the data management back-end will need additional functionality as the dynamic, on-demand approach to the production of GAR and country reports will be adapted. In the absence of a clear two-year cycle of report production, issues such as persistent identifiers, versioning of the data and metadata need to be reviewed. Initial assessment of solutions such as KIT Datamanger [START_REF] Wm | Universität des Landes Baden-Württemberg and nationales Forschungszentrum in der Helmholtz-Gemeinschaft: Kit Data Manager -The Research Data Repository Platform[END_REF] and CKAN [START_REF]CKAN -The open source data portal software[END_REF] as tools to meet these new requirements is already ongoing.

In response to the computing challenges related to risk calculation outlined in the chapter 4.1, more advanced parallelisation approaches are already being studied by the collaboration. For example, both MapReduce [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF] and MPI [START_REF][END_REF] based approaches could allow parallelisation beyond the limitations of the shared memory space andwith sufficient resources -even allow in-memory processing of the data.

The new data management approaches and tools discussed will obviously also have an impact on the processes that are internal to GAR collaboration manages the data. They may also influence the interfaces the GAR data services can support for the thirdparty analysis. As an example of a potential explorative research topic, investigating approaches where the analysis of the data could be performed at the storage location in a flexible manner (e.g. by shipping the analysis code embedded in a virtual machine or software container to data) could support more efficient and flexible management and use of growing and increasingly dynamic GAR data.

Fig. 1 .

 1 Fig.1. A small sample of some of the risk results of Germany as presented in the online version of the GAR. The average annual losses are put in the context of key macroeconomic indicators of the country, such as the Gross Fixed Capital Formation (GFCF) that measures the annual net increase of fixed assets owned by business sector, government and households and Social Expenditure. This contextualisation illustrates the load on the society natural disasters represent more effectively than mere average annual loss figures would.

Fig. 2 .

 2 Fig. 2. Part of the Annual Average Loss data of Guatemala illustrate how the impact on society is considerably more severe despite the absolute losses being lower than the German ones (Figure 1)

Fig. 5 .

 5 Fig. 5. A screenshot of an early version of the download portal (proof of concept) used to gather feedback related to the functionality.


Sacha Gobeyn
email: sacha.gobeyn@ugent.be

Peter L M Goethals

A variable length chromosome genetic algorithm approach to identify species distribution models useful for freshwater ecosystem management

Keywords: species distribution models, model identification, genetic algorithms, freshwater management, macroinvertebrate species, input variable selection, parameter estimation

Increasing pressure on freshwater ecosystems requires river managers and policy makers to take actions to protect ecosystem health. Species distribution models (SDMs) are identified as appropriate tools to assess the effect of pressures on ecosystems. A number of methods are available to model species distributions, however, it remains a challenge to identify well-performing models from a large set of candidate models. Metaheuristic search algorithms can aid to identify appropriate models by scanning possible combinations of explanatory model variables, model parameters and interaction functions. This large search space can be efficiently scanned with simple genetic algorithms (SGAs). In this paper, we test the potential of a variable length chromosome SGA to perform parameter estimation (PE) and input variable selection (IVS) for a macroinvertebrate SDM. We show that the SGA is an appropriate tool to identify fair to satisfying performing SDMs. In addition, we show that SGA performance and the uncertainty varies as a function of the chosen hyper parameters. The results can aid to further optimise the algorithm so models explaining species distributions can be identified and used for analysis in river management.

Introduction

Freshwater ecologist and river managers are in need for system analysis techniques to investigate a wide range of ecological questions and support decision making. Species distribution models (SDMs) aiming to describe the species response to driving processes, have shown to be valuable tools in ecosystem health management. Many approaches to identify SDMs are available, however, the challenge remains to test a large set of candidate explanatory models. Genetic algorithms (GAs) classified under evolutionary algorithms and inspired by various mechanisms observed in evolution (i.e. reproduction, mutation, selection) are promising approaches to evaluate a large search space [START_REF] Gibbs | A genetic algorithm calibration method based on convergence due to genetic drift[END_REF][START_REF] Hamblin | On the practical usage of genetic algorithms in ecology and evolution[END_REF][START_REF] Maier | Evolutionary algorithms and other metaheuristics in water resources: Current status, research challenges and future directions[END_REF]. Consequently, GAs are used to select input variables (input variable selection, IVS) for SDMs by using them as a wrapper for data-driven approach [START_REF] D'heygere | Use of genetic algorithms to select input variables in decision tree models for the prediction of benthic macroinvertebrates[END_REF]. They are also used to estimate parameter values (parameter estimation, PE) for fuzzy logic SDMs [START_REF] Fukuda | Effect of model formulation on the optimization of a genetic Takagi-Sugeno fuzzy system for fish habitat suitability evaluation[END_REF][START_REF] Van Broekhoven | Interpretability-preserving genetic optimization of linguistic terms in fuzzy models for fuzzy ordered classification: an ecological case study[END_REF]. PE and IVS are important aspects in SDM identification and it can hypothesized whether a joint approach can be encoded in GAs.

In this paper, we present the use of a simple genetic algorithm (SGA) for PE and IVS for an SDM. To do so, we encode the optimisation problem in a variable length chromosome. The approach is tested for a freshwater species, cloeon dipterum, with the Limnodata of the Netherlands. The acquired SDM performance, parameters and input variables are analysed. In addition, a sensitivity analysis is done to test the effect of the algorithm hyper parameters on the SGA performance (Section 3). The results of this approach are discussed in section 4.

SDM development

The SDM is developed by following a four step approach (Figure 1). First, a number of ecological concepts are used to define the model. Second, the data are gathered and processed to construct the model (step 3). In a final step, a search algorithms is implemented and used to identify well-performing models.

Model concept

Filter theory is used as basic concept for the SDM. In this theory, the realized species assemblage is explained by a number of hierarchical filters, i.e. dispersal, abiotic and biotic filters [START_REF] Guisan | SESAMa new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages[END_REF][START_REF] Poff | Landscape filters and species traits: towards mechanistic understanding and prediction in stream ecology[END_REF]. Here, it is used because of its structural nature dividing the explanatory processes of species presence/absence in several filters. Only abiotic filtering is considered because the effect of pollutants on the species assemblage is assumed as the most relevant source of information for ecosystem health management. Species response curves (or habitat preference curves) defining the biological response to abiotic gradients are used as to reflect the abiotic filters. The biological response can be expressed by many measures, i.e. species presence, abundance, density, usable area or volume. Species presence is used as a measure for biological response because it is assumed to be a robust measure for biological response [START_REF] Fukuda | Abundance versus presence/absence data for modelling fish habitat preference with a genetic Takagi-Sugeno fuzzy system[END_REF].

In this paper, fine-scale and large-scale abiotic filters are considered. Fine-scale filters are filters acting at a local scale filtering species due to point specific pollution. In addition, the river typology characteristics (e.g. geology, river/catchment slope, ..) are considered to be large-scale filters which act on a river or catchment scale [START_REF] Poff | Landscape filters and species traits: towards mechanistic understanding and prediction in stream ecology[END_REF].

Data

The Limnology Neerlandica database (http://www.stowa.nl/) and information on the river typology [START_REF] Elbersen | Definitiestudie kaderrichtlijn water (KRW) I. Typologie nederlandse oppervlaktewateren[END_REF] are processed and compiled to a coupled database. The Limnodata is a database containing observations of the biology (macroinvertebrates, fish and macroflora) and physical-chemical state over 20 years in the Netherlands. The river typology is defined as a function of river catchment characteristics, i.e. average river slope, water source, average river width, catchment area, tidal influence, catchment geology [START_REF] Sandin | Stream and river typologies -major results and conclusions from the STAR project[END_REF].

The observations of the macroinvertebrate species cloeon dipterum are extracted from the Limnodata. The records are transformed from abundance to presence/absence in order to get an insight in the spatial and temporal distribution patterns of the species. Outliers in the physico-chemical data are investigated by inspecting summary statistics (mean, minimum, maximum and percentile values) and visually analysing box plots, histograms and dot plots. A number of variables are tested to physical boundaries. For instance, the width and depth of rivers are assessed as a function of the river type. In addition, the mass balance for nitrogen and phosphorus is inspected. In total 133 values are inspected in-depth leading to the omission of 102 records from the data. Finally, the correlation between variables is calculated so to exclude highly correlated variables and reduce dimensionality of the problem (Table 1).

Model construction

Species response curves are defined for the fine-scale filters (continuous variables). The species response curves are assumed to have a non-symmetric unimodal trapezoid shape chosen as a simplification of a bell-shaped curve [START_REF] Austin | Species distribution models and ecological theory: a critical assessment and some possible new approaches[END_REF][START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF]. The curves are allowed to be asymmetric so they can skew from extreme (heavy polluted) conditions [START_REF] Austin | Species distribution models and ecological theory: a critical assessment and some possible new approaches[END_REF][START_REF] Hirzel | Habitat suitability modelling and niche theory[END_REF]. Four parameters (𝑎 1 , 𝑎 2 , 𝑎 3 and 𝑎 4 ) are used to define the trapezoid curve:

𝑆𝐼 𝑓 (𝑥 𝑖 𝑗 ) = { 0 if 𝑥 𝑖 𝑗 < 𝑎 1 (𝑥 𝑖 𝑗 -𝑎 1 ) (𝑎 2 -𝑎 1 ) if 𝑥 𝑖 𝑗 ∈ [𝑎 1 , 𝑎 2 [ 1 if 𝑥 𝑖 𝑗 ∈ [𝑎 2 , 𝑎 3 ] (𝑎 4 -𝑥 𝑖 𝑗 ) (𝑎 4 -𝑎 3 ) if 𝑥 𝑖 𝑗 ∈ ]𝑎 3 , 𝑎 4 ] 0 if 𝑎 4 < 𝑥 𝑖 𝑗 (1)
With 𝑆𝐼 𝑓 , the suitability index for the fine-scale filters, 𝑥 𝑖 𝑗 , the input value 𝑖 (∈ {0,1, . . , 𝑁}, 𝑛 data points) for variable 𝑗 (∈ {0,1, . . , 𝑀}, 𝑚 variables). The parameters 𝑎 1 and 𝑎 4 describe the range of the conditions in which a species is able to survive. The parameters 𝑎 2 and 𝑎 3 describe the preferable range of conditions for the species (i.e. 𝑆𝐼 =1). The values of 𝑎 1 and 𝑎 4 are set by the minimum and maximum values of the observations for which the species is observed. For the large-scale abiotic filters, suitability indices are defined based on a set of parameters (𝑎 1 , 𝑎 2 , .., and 𝑎 r ) and the class (categorical):

𝑆𝐼 𝑙 (𝑥 𝑖 𝑘 ) = { 𝑎 1 if 𝑥 𝑖 𝑘 = C 1 𝑘 𝑎 2 if 𝑥 𝑖 𝑘 = C 2 𝑘 … 𝑎 r if 𝑥 𝑖 𝑘 = C 𝑟 𝑘 (2)
With 𝑆𝐼 𝑙 , the suitability index for the large-scale filter, 𝑥 𝑖 𝑘 , the input value 𝑖 (∈ {0,1, . . , 𝑁}, 𝑛 data points) for categorical variable 𝑘 (∈ {0,1, . . , 𝑂}, 𝑜 variables). The habitat suitability index (𝐻𝑆𝐼) for a point 𝑖 is calculated by multiplying the geometric mean for the fine and large-scale filters:

𝐻𝑆𝐼 𝑖 = (∏ 𝑆𝐼(𝑥 𝑖 𝑗 ) 𝑚 𝑗=1 ) 1 𝑚 ⁄ * (∏ 𝑆𝐼(𝑥 𝑖 𝑘 ) 𝑜 𝑘=1 ) 1 𝑜 ⁄ (3)

Model identification with simple genetic algorithms

The aim of the model identification tool is to identify a number of input variables and coupled species response curve parameters with an optimisation algorithm. This algorithm has to be able to efficiently search a large unconstrained space since it is difficult to a priori define the shape of a species response (skewed, Gaussian, …). In addition, a number of solutions is possibly more informative than one solution. Therefore, it is preferred to obtain an ensemble. The tool should be compatible with high performance computing to facilitate repeated runs for uncertainty analysis. Even more, it is required to be an open source package, available freely online, so to increase code and approach transparency (Figure 2). An SGA with three operators, i.e. selection, crossover and mutation is implemented and used as optimisation algorithm. An SGA requires the encoding of the phenotype, i.e. the model, in a genotype. This genotype is typically coded as a binary string. This string is then translated to a model in a genotype-phenotype mapper. A list of lists is programmed to implement a variable length chromosome (Figure 3). The genome is defined by a second order binary string when a bit in the first order binary string has a value of one. The first order binary string is translated in a mapper by either in-or excluding the variable (one = present, zero =absent). The second order binary string is translated to parameter values of 𝑎 2 and 𝑎 3 in the mapper function by transforming every three bit sequence to an integer representation which is used to define the values for 𝑎 2 and 𝑎 3 (equation 1) and 𝑎 r (equation 2). The tournament selection method is used to select the fittest individuals from a population as parents [START_REF] Goldberg | A comparative analysis of selection schemes used in genetic algorithms[END_REF]. The selection rate defined as the fraction of the population that survives for the next step of mating is multiplied with the population size to obtain a number of parents. In the crossover operator, the parents are randomly paired to mate and produce offspring with a certain rate, i.e. crossover rate. If mating does not occur, the parents are replaced in the population. The last operator, mutation, is defined as the probability that a random gene is assigned a new value (0 → 1 or 1 → 0). The fitness of the chromosomes is the sum of squared errors (SSE) calculated with the 𝐻𝑆𝐼 𝑖 values and the observed presence or absence 𝑃𝑟 𝑖 :

𝑆𝑆𝐸 = ∑ ( 𝑛 𝑖=1 𝑃𝑟 𝑖 -𝐻𝑆𝐼 𝑖 ) 2 (4)
3

Results

The SGA is implemented and used to identify near-optimal models for the species cloeon dipterum. In the first part of this section, a set of hyper parameters (mutation and crossover rate) for the SGA are tested so to estimate the effect of hyper parameter choice on the algorithm performance. In the second part, the results found with the SGA and near-optimal hyper parameters are used to analyse the acquired model structure and performance.

Sensitivity of SGA

The SGA sensitivity as a function of the hyper parameter values are shown in Figure 4.

For this experiment, an initial near-optimal set of parameters is determined by following the guidelines of [START_REF] Gibbs | A genetic algorithm calibration method based on convergence due to genetic drift[END_REF]. The required number of chromosomes 𝑃 are estimated by applying equation 5:

𝐹𝐸 𝑃 log 10 (1 - 1 𝑃 ) = -𝑀 -log 10 (√ 𝑙 12 ) (5)
With 𝑀 equal to three, 𝐹𝐸, the number function evaluations determined by dividing the computational time available by the average runtime of one simulation and 𝑙, the chromosome length. For 𝑙, the maximum possible length of the chromosome is used ( = 111 = three bits * (two parameters * 12 continues variables + 13 parameters for categorical variables)). With equation 5, 100 is found as a value for 𝑃. The mutation rate is calculated by dividing five by 𝑃 (pm = 0.05*100 %) and the crossover rate (pc) is set to 100 % [START_REF] Gibbs | A genetic algorithm calibration method based on convergence due to genetic drift[END_REF]. It is assumed that the performance of the SGA is near-optimal with these values. In order to verify the choice of the values, the sensitivity of the SGA performance to the values is checked by assessing the effect of the surrounding values of the found near-optimal values for the crossover and mutation rate (nine point grid with pc = {60, 85, 100} and pm = {1, 5, 10}).

The best found solution follows a similar evolution for the nine sets of hyper parameters with a varying convergence and performance. When inspecting the effect of the mutation rate (pm) on the performance of the algorithm, one observes that the SGA analysis with a mutation rate of 1 % gives on average the best solutions (Figure 4, left panel). The evolution of the best solution found with a mutation rate of 5 % is similar whereas a higher mutation rate (10 %) leads to less optimal solutions. The initial speed with which these solutions are found is highest for a mutation rate of 5 %, however, the population converges -on average -earlier. For the crossover rate (Figure 4, right panel), one observes that the sensitivity of the performance is lower than for the mutation rate. One observes that a varying degree of uncertainty is observed for different hyper parameter values (Figure 4 and5). The uncertainty is estimated by repeatedly running the SGA with a number of initial conditions and preserving the best solution over the generations for every SGA run. The variation of this uncertainty follows a hyperbole as a function of the mean SSE (and thus the generation) (Figure 5, left panel). At the point of convergence (low SSE, Figure 5, right panel), the uncertainty on the found nearoptimal solutions for a crossover rate of 85 % is lower than for a crossover rate of 100 %. This seems to suggest that the crossover rate of 85 % is an appropriate choice to reduce SGA analysis uncertainty.

Analysis of identified SDMs

The acquired models with the SGA (pm= 1 %, pc = 85 %) are evaluated by calculating the Cohen's Kappa (Kappa) and area under the receiver operator curve (AUC) (see [START_REF] Mouton | Ecological relevance of performance criteria for species distribution models[END_REF] for mathematical description). The acquired models are assessed to have a fair to satisfying performance. The mean Kappa is equal to 0.33±0.03 which is assessed as fair (Kappa ∈ [0.2,0.4], see [START_REF] Gabriels | Analysis of macrobenthic communities in Flanders, Belgium, using a stepwise input variable selection procedure with artificial neural networks[END_REF]). The mean AUC is equal to 0.7±0.03 which is assessed as satisfying (AUC > 0.7, see [START_REF] Manel | Evaluating presence-absence models in ecology: the need to account for prevalence[END_REF]). In Figure 6, the model structure and accompanied uncertainty found by repeatedly running the SGA is shown. The support (%) for a model variable is calculated as a measure of variable importance by dividing the number of times a variable is selected by the SGA by the total number of SGA analysis. The support for the variable pH is very high (99 %) whereas the support for the river slope, catchment area, tides and geology is lower and uncertain. In Figure 7, the species response curves and the accompanied uncertainty for the variable pH is shown. Either a response with very steep boundaries or a triangular response is observed. The uncertainty is shown for three values of the mutation rate (constant crossover rate). It is observed that the uncertainty on the acquired curves increases for higher mutation rate. This patterns is similar to the increase of uncertainty in the convergence of the objective function (Figure 5, right panel). When inspecting the uncertainty on the parameters of the categorical variables (not shown here), one observes a rather high uncertainty. In conclusion, the uncertainty in the objective function is reflected in the uncertainty of the model structure. Fig. 6. Support for variable inclusion for repeated SGA analysis (pc = 85 %, pm = 1 %). The support is calculated by dividing the number of times a variable is selected by the SGA by the total number of analysis (i.e. 100). The uncertainty is estimated with the Shannon entropy [START_REF] Sandin | Stream and river typologies -major results and conclusions from the STAR project[END_REF].

Discussion and outlook

In this paper, a variable length chromosome SGA is implemented and used to jointly perform IVS and PE. The implemented algorithm is able to identify fair to satisfying models. The uncertainty on the acquired species response curve parameters is rather low, at least for the variable with a high support. In addition, it is observed that the uncertainty on the acquired near-optimal solution is not equal over different values of the mutation and crossover rate. The accuracy of the models could be improved by increasing the precision of the binary encoding used for the algorithm. In the current implementation every three bits code one parameter of the species response curves (see Figure 3). This allows to encode eight discrete values for every parameter. The representation restricts the possible parameter values to a limited set defined by the lower and upper boundary of the parameter interval and the number of bits [START_REF] Van Broekhoven | Interpretability-preserving genetic optimization of linguistic terms in fuzzy models for fuzzy ordered classification: an ecological case study[END_REF]. Increasing the number of bits for the binary encoding might increase the precision but will also increase the length of the chromosome. Consequently, different near-optimal values for the hyper parameters will be obtained with equation 5. When testing the required number of chromosomes, for a fixed number of 𝐹𝐸, one observes that the found number of chromosomes (and thus mutation rate, see [START_REF] Gibbs | A genetic algorithm calibration method based on convergence due to genetic drift[END_REF]) is almost equal for higher chromosome lengths. For example, for a bit length of three, a maximum chromosomes length of 111 (three bits * (two parameters * 12 continues variables + 13 parameter for categorical variables)) leads to a population size of 112, whereas for a four and six bit problem (𝑙 = 148, 𝑙 = 185) a number of 111 and 110 chromosomes is found. Since the determined near-optimal values for the hyper parameters for varying chromosome lengths does not vary, it is expected that the performance and uncertainty of the SGA will not vary as a function of the length used to encode the optimisation. This suggests that increasing the precision of the binary encoding will not influence the performance and uncertainty of the SGA analysis. Additional experiments with the SGA should confirm this hypothesis.

A hyperbolic relation is found between the uncertainty on the SGA analysis and the found near-optimal solution. At the start of the analysis, the uncertainty is rather small, and increases with the number of generations to finally converge to a value as the SGA converges. There are differences in the amount of uncertainty at convergence for varying values of mutation and crossover rate. For low mutation rates, the uncertainty on the found near-optimal solution declines as the crossover rate is lowered. For higher mutation rates, this relation is inverse but less apparent. In general the guidelines by [START_REF] Gibbs | A genetic algorithm calibration method based on convergence due to genetic drift[END_REF] are assessed as appropriate for these type of problems, since with these settings the SGA is able to reduce the prediction error of the models (mean SSE declines from approximately 140 to 40). Options to further improve the algorithm performance can be to improve the exploitive character of the algorithm by combing the genetic algorithm with a hill climbing (HC) approach or to vary the mutation and crossover rates over the generations. Further research can investigate whether these implementation have a significant added value for SDM identification and whether they can reduce the uncertainty of the analysis.

Genetic algorithms have shown to be valuable for PE and IVS in species distribution modelling [START_REF] D'heygere | Use of genetic algorithms to select input variables in decision tree models for the prediction of benthic macroinvertebrates[END_REF][START_REF] Fukuda | Effect of model formulation on the optimization of a genetic Takagi-Sugeno fuzzy system for fish habitat suitability evaluation[END_REF]. In this study, a variable length chromosome implementation of an SGA is presented to jointly perform PE and IVS. The results tested for one species are promising, however, it should be further investigated how the performance of the algorithm varies as function of the algorithm settings. In addition, the approach should be validated by applying the SGA for different species.

The current available software is an open source package implemented in the Python programming language [START_REF] Gobeyn | Species distribution model identification tool (SDMIT)[END_REF][START_REF] Gobeyn | Input variable selection with a simple genetic algorithm for conceptual species distribution models: A case study of river pollution in Ecuador[END_REF]. Many other packages are available (Generalized Linear Models, GLM, in the R programming language or Genetic Algorithm for Rule set Production/Prediction (GARP) software [START_REF] Stockwell | The GARP modelling system: problems and solutions to automated spatial prediction[END_REF]). For instance, the GLM R package is an userfriendly package useful for ecologist, however, automated running a number of analysis to estimate uncertainty is difficult. Even more, the statistical approaches present a number of boundary conditions to the shape of the species response. In the developed approach, these boundary conditions are relaxed as is the case for GARP. The difference with GARP is that the SDMIT approach is designed to run on high performance clusters whereas GARP was initially designed for single-run analysis in a graphical user interface environment. In addition, machine learning approaches like decision trees and support vector machines are available [START_REF] Fukuda | Habitat prediction and knowledge extraction for spawning European grayling (Thymallus thymallus L.) using a broad range of species distribution models[END_REF]), however, the disadvantage of these approaches is that the tools are not implemented specifically for the optimisation of SDMs and thus often lack the ecological theoretical background. Consequently they are used as data mining approaches rather than model optimisation algorithms. The SDMIT packages is an answer to these limitations. With this, SDMs can obtained that improve the insight in species and community response to environmental changes.

Fig. 1 .

 1 Fig. 1. Overview of methodology to develop SDMs (adapted after [2]).

Fig. 2 .

 2 Fig. 2. Requirements for model identification tool for species distribution models.

Fig. 3 .

 3 Fig.3. Definition of chromosome and mapper function. The genome is programmed as a list of lists, where a second order binary string is defined when a bit of the first order string has the value of one. Every three bits of this second order string are translated to an integer which is used to define the values of the parameters 𝑎 2 , 𝑎 3 and 𝑎 r . In this example, the second and sixth variable are considered in the model. The parameters for the species response curves are defined by second order binary strings (six bits). The first three bits for variable X2 are used to define aX2,2 and the last three bits to define aX2,3. A binary coding is used to define a fraction (i.e. 5/7 and 6/7) of the total range r (a4 -a1) which is added to the parameter a1 to obtain values for a2 and a3. Parameters a2 and a3 are respectively bounded by the range [a1, a3] and [a2,a4]. For the categorical variables, the parameters 𝑎 𝑟 are bounded by zero and one.

Fig. 4 .

 4 Fig. 4. Evolution of SSE (inverse of fitness) as a function of the number of generations. On the left, the results are shown for varying mutation rates and constant crossover rates. On the right, the results are shown for constant mutation rates and varying crossover rates. The uncertainty on the analysis is acquired by repeating the SGA a number of times with different initial conditions and preserving the best solution every generation.

Fig. 5 .

 5 Fig. 5. Standard deviation on SSE as a function of the mean value of the SSE, for varying values of the hyper parameters (%). The right panel zooms in a narrower range of the left panel.

Fig. 7 .

 7 Fig. 7. Uncertainty on acquired species response curves with the repeated SGA analysis. From left to right, a mutation rate of 1, 5 and 10 % is used (constant crossover rate = 85 %).

Table 1 .

 1 Overview

		of physico-chemical variables. #n: not included because of insufficient sam-
	ples after coupling with biological data. ex. = excluded, corr. = correlated to, r = spearman rank
	correlation, min. = minimum, 𝑋 ̃ = median, 𝑋 ̅ = mean, max. = maximum, Chlor. a = chlorophyll
	a, Cond. = conductivity, Transp. = transparency, Kjel. = kjeldahl N, R. = river, Temp. = temper-
	ature.
	variable	ex.	reason	min.	𝑿 ̃	𝑿 ̅	max.
	%DO	x	Corr. DO	0.00	80.00 78.07	277.00
			(r = 0.89)
	BOD2	x	#n	10.00	10.00 92.87	2000.00
	BOD5			0.05	2.00	3.55	360.00
	Chloride			1.00	40.00 56.99	1250.00
	Chlor. a	x	#n	0.10	9.00	19.51	1170.00
	COD	x	#n	2.00	26.00 32.57	200.00
	Cond.	x	Corr. Chloride	0.50	50.00 52.48	542.00
			(r = 0.79)
	DO			0.00	8.80	8.64	29.00
	Transp.			0.00	0.50	0.50	3.00
	Flow	x	#n	0.00	0.10	0.69	33.34
	Kjel. N	x	Corr. NH4	0.00	1.70	2.60	70.00
			(r = 0.91)
	NH3-N			0.00	0.01	0.03	6.10
	NH4-N			0.00	0.40	1.13	80.00
	NO2-N			0.00	0.06	0.10	6.30
	NO3-N			0.00	3.50	4.89	64.00
	PO4-P			0.00	0.07	0.28	26.00
	pH			3.60	7.40	7.33	10.40
	R. depth	x	#n	0.00	0.40	0.66	5.00
	R. width	x	#n	0.02	3.00	5.94	135.00
	SO4			1.00	62.00 68.15	6200.00
	Temp.			-1.00	11.50 11.73	32.00
	Total N	x	Corr. NO3-N	0.05	5.56	7.07	66.30
			(r = 0.93)
	Total P	x	Corr. to PO4-P	0.00	0.20	0.47	29.00
			(r = 0.92)
	Velocity	x	#n	0.00	20.00 24.15	300.00

Acknowledgments. Sacha Gobeyn is supported by a Bijzonder Onderzoeksfonds (BOF) project related to the Ecuador Biodiversity Network of the Vlaamse Interuniversitaire Raad-Universitaire Ontwikkelingssamenwerking (VLIR-UOS). This research was performed in the context of the VLIR Ecuador Biodiversity Network project. The computational resources (Stevin Supercomputer Infrastructure) and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by Ghent University, the Hercules Foundation and the Flemish Government -department EW.