SuperflexPy 1.2.0: an open source Python framework for building, testing and improving conceptual hydrological models

10 Catchment-scale hydrological models are widely used to represent and improve our understanding of hydrological processes, and to support operational water resources management. Conceptual models, where catchment dynamics are approximated using relatively simple storage and routing elements, offer an attractive compromise in terms of predictive accuracy, computational demands and amenability to interpretation. This paper introduces SuperflexPy, an open-source Python framework implementing the 15 SUPERFLEX principles (Fenicia et al., 2011) for building conceptual hydrological models from generic components, with a high degree of control over all aspects of model specification. SuperflexPy can be used to build models of a wide range of spatial complexity, ranging from simple lumped models (e.g. a reservoir) to spatially distributed configurations (e.g. nested sub-catchments), with the ability to customize all individual model elements. SuperflexPy is a Python package, enabling modelers to exploit 20 the full potential of the framework without the need for separate software installations, and making it easier to use and interface with existing Python code for model deployment. This paper presents the general architecture of SuperflexPy, and illustrates its usage to build conceptual models of varying degrees of complexity. The illustration includes the usage of existing SuperflexPy model elements, as well as their extension to implement new functionality. SuperflexPy is available as open-source code, and 25 can be used by the hydrological community to investigate improved process representations, for model comparison, and for operational work. A comprehensive documentation is available online and provided as supplementary material to this paper. https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c © Author(s) 2020. CC BY 4.0 License.


Conceptual hydrological models
Catchment-scale hydrological models are widely used to predict catchment behavior under natural and human-impacted conditions, as well as to represent and improve our understanding of internal catchment 70 functioning (e.g. Beven, 1989). For example, catchment models underlie projections of climate change impact on groundwater recharge and streamflow (e.g., Eckhardt and Ulbrich, 2003), are used as tools for hypothesis testing to identify dominant hydrological processes (e.g., Clark et al., 2011b;Hrachowitz et al., 2014;Wrede et al., 2015), and are used to inform agricultural practices such as irrigation scheduling (e.g., McInerney et al., 2018) and pesticide application (e.g., Moser et al., 2018;Ammann et al., 2020). The 75 typical use of hydrological models is to simulate or forecast the streamflow response (runoff) of a catchment to rainfall; for this reason they are often referred to as rainfall-runoff (RR) models (e.g., Moradkhani and Sorooshian, 2009). However, their application extends to the simulation of other environmental variables, including internal catchment variables such as groundwater levels (e.g., Seibert and McDonnell, 2002) and soil moisture (e.g., Matgen et al., 2012), as well as water chemistry (e.g., 80 Bertuzzo et al., 2013;Ammann et al., 2020) An important class of catchment models are "process based" models, which describe the cascade of processes transforming catchment inputs (e.g. precipitation) into outputs (e.g. streamflow) though conservation equations. These models are an appealing choice due to their broad physical underpinnings, as well as their ability to represent internal catchment processes and potential for predicting catchment 85 responses under changing environmental conditions. Process based models can be classified according to the nature of their constitutive equations (e.g. conceptual or physically based) and their spatial resolution (e.g. lumped or distributed) (e.g., Refsgaard, 1996).
In terms of spatial resolution, conceptual models can be applied in a lumped configuration (treating the entire catchment as a single unit) if the interest is in modeling integrated catchment outputs (e.g. 95 streamflow). Alternatively, distributed configurations can be used when the interest is in modeling hydrological behavior at individual landscape sections (e.g., sub-catchments). In such distributed setups, the catchment is subdivided into spatial elements such as sub-catchments (e.g., Feyen et al., 2008;Lerat et al., 2012), Hydrological Response Units (HRUs) (e.g., Arnold et al., 1998;Fenicia et al., 2016;Dal https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License. Molin et al., 2020b), or grids (e.g., Samaniego et al., 2010). A common strategy for developing distributed 100 conceptual models is to represent individual landscape elements using independent (non-interacting) lumped models, and then obtain total catchment outflow by aggregating the outflows from these individual models, potentially incorporating flow routing elements to represent routing delays. This strategy is often referred to as "semi-distributed" modelling, and typically employs discretization based on principles of "hydrological similarity" (e.g., Sivapalan et al., 1987), such as HRUs (e.g., Leavesley, 105 1984). In many cases, semi-distributed modelling achieves good predictive ability while greatly simplifying model representation and reducing computational demands compared to fully-integrated 2D/3D distributed models such as Parflow (Maxwell, 2013) or Mike She (Refsgaard and Storm, 1995), which typically use much smaller landscape elements and explicitly model lateral exchanges. For the purposes of this presentation, we consider semi-distributed modelling to be a special case of distributed 110 modelling.

Hydrological model structure and flexible modeling frameworks
The selection of model structure has preoccupied researchers and practitioners since the early days of hydrological modelling (e.g., Ibbitt and O'Donnell, 1971;Moore and Clarke, 1981;Jakeman and Hornberger, 1993). Although in principle the physical laws governing hydrological processes are the 115 same everywhere, the diversity of catchment conditions in terms of topography, soil, geology, vegetation, and anthropogenic influence, results in remarkably different manifestations of these physical laws at the catchment scale. These local differences, also termed "uniqueness of place" (Beven, 2000), considerably limit our ability to develop generalizable hydrological hypotheses (e.g., Wagener et al., 2007).
Model structure selection has led to multiple research directions, including the search for a single model 120 structure that achieves good prediction across all catchments (the "fixed" model paradigm), and the search for model structures best suited for particular types of environmental conditions (the "flexible" model paradigm). Whether in search of a single model or multiple models, model selection necessarily relies on a process of model development, comparison, and refinement. Approaches to formalize this process include the top-down approach (e.g. Sivapalan et al., 2003), the system identification approach (e.g 125 Young, 1998), and the method of multiple working hypotheses (e.g., Clark et al., 2011a). These approaches are not mutually exclusive, as the idea of comparing different model representations is ubiquitous in model development and empirical science in general.
The process of model development, comparison, and refinement can be facilitated using flexible modeling frameworks, which enable hydrologists to hypothesize, implement, and (eventually) test and refine 130 different model structures. Flexible frameworks have themselves developed along multiple directions.
In this paper, we focus on flexible frameworks intended for conceptual hydrological modeling. Many frameworks have been developed for such purpose, and offer different degrees of flexibility. For example, FUSE (Clark et al., 2008) allows exchanging the components of 4 common models (Sacramento, PRMS, TOPMODEL and ARNO/VIC). SUPERFLEX  allows building model structures 140 from generic elements (reservoirs, lag function and connections). CMF (Kraft et al., 2011), represents the model as an abstract network of elements and can be adopted for conceptual models. PERSiST (Futter et al., 2014) allows to create semi-distributed bucket-type models and is designed to be coupled with a water quality model. ECHSE (Kneis, 2015) is a framework for development of object-based conceptual hydrological model engines. MARRMoT (Knoben et al., 2019) provides a unified implementation of 46 145 existing conceptual models (including GR4J, HBV, and others).
When discussing any mathematical model, it is relevant to distinguish its conceptual principles from its software implementation. For example, common conceptual models, such as GR4J, HBV and TOPMODEL, exist in many public and in-house versions and in many computer languages (Excel, R, Matlab, Fortran, C, and others). FUSE, to our knowledge, has implementations in Fortran (Clark et al., 150 2008) and R (Vitolo et al., 2016). SUPERFLEX, besides its original Fortran implementation , has also been coded in Matlab (David et al., 2019).
Ideally, the software implementation of a flexible framework should fulfill its goals, i.e., (1) cover the envisioned range of applications (e.g., flexibility, spatial extension, processes representation, etc.), and (2) achieve this in a sufficiently "practical" way (i.e., it should not require complicated and abstract setup 155 procedures to define its configuration).
For conceptual models, a flexible model framework should arguably cover the following "realms": 1. Lumped models; 2. Distributed setups, including simulation of sub-catchments and flows/processes at internal points; 3. Substance transport modelling, including water isotopes, pesticides, etc; 160 4. Ability to reproduce existing models, when necessary.
1. Ease of use, including installation, learning, and operation. Interoperability with external software, for example for model calibration and uncertainty analysis, is of obvious relevance because hydrological models are often used as parts of larger-scale projects and operations. 165 2. Ease of modifications and extensions. Even a comprehensive software implementation will eventually require extension. For example, a modeling framework intended for streamflow simulation may require extension to simulate water chemistry. Another type of modification might be a switch to a different numerical implementation better suited for parallel computing, etc.
3. Computational efficiency. Hydrological model applications, especially including calibration and 170 uncertainty quantification, may require thousands or even millions of model runs.
Arguably, these requirements are not collectively fulfilled by available software implementations of flexible frameworks. For example, only CMF and ECHSE provide full flexibility in terms of model structure selection. In some frameworks, the intended flexibility is obtained by enabling and tuning the components of a master structure (e.g. FUSE, PERSiST). The original implementation of SUPERFLEX 175 in Fortran , used for research case studies though not released publically, also used this "master structure" approach. In other flexible frameworks, the user can chose from an extensive library of existing model configuration (e.g. MARMMoT). Some flexible frameworks are limited to lumped configurations (e.g. FUSE, MARMMoT). Substance transport is currently partially possible with PERSiST, by interfacing with additional software. 180 Clearly, achieving all these objectives simultaneously is challenging, and entails juggling multiple obvious and less obvious tradeoffs. For example, the intended flexibility of a framework may come at the expense of ease of use, similar to how computer languages have varying degrees of abstraction from the hardware behavior. Implementing a practical flexible framework therefore requires careful code design, experimentation, and inevitably, some compromises. 185 This work pursues the flexible framework objectives by building upon the concept of SUPERFLEX Kavetski and Fenicia, 2011;Fenicia et al., 2014;Fenicia et al., 2016). A key attractive feature of SUPERFLEX as a modelling concept is the fine "granularity" of model structures it can support, which enables systematic and detailed hypothesis testing . For example, the hydrologist should have the ability to change the functional form of a single flux expression in one of the 190 model elements, as well as to change such flux expressions in multiple specific parts of the model.
The development of the proposed framework capitalizes on the authors' collective experience in using the earlier implementations of SUPERFLEX in a series of empirical case studies over the last decade, ranging from lumped model implementations (e.g., Kavetski and Fenicia, 2011;Fenicia et al., 2014), to distributed setups (e.g. Fenicia et al., 2016;Dal Molin et al., 2020b), interpretation in the context of 195 https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License. fieldwork insights (e.g., Wrede et al., 2015), large scale model intercomparisons (e.g., van Esse et al., 2013), and the inclusion of pesticide/substance transport (e.g. Ammann et al., 2020). These applications have highlighted the versatility of SUPERFLEX principles to solve increasingly complex modelling problems, and have led to insights into software design and configuration aspects not available in the earlier implementations. This study reports on these developments and offers an open source 200 implementation of SUPERFLEX for use by the hydrological community.

Aims
This work introduces SuperflexPy, an open-source Python software that implements the principles of the SUPERFLEX framework for conceptual hydrological model development. Our objectives are as follows:

Present SuperflexPy and its basic building blocks (components): elements, units, nodes, and 205
network.
2. Illustrate how SuperflexPy can help hydrologists implement a conceptual model structure at the desired level of internal complexity and spatial resolution -including recreating existing models or developing new ones.
3. Provide a broad discussion of how the SuperflexPy contributes to the toolkits available to the 210 hydrological community, including existing flexible frameworks, in terms of intended scope of application, advantages, and limitations.
The paper is organized as follows: Sect. 2 presents SuperflexPy to the hydrological community; Sect. 3 illustrates selected applications of the framework including the setup of SUPERFLEX configurations used in earlier case studies, as well as how to use SuperflexPy to create new elements; Sect. 4 provides 215 SuperflexPy implementation details useful for understanding the usage and general potential of the framework; Sect. 5 discusses the scope of SuperflexPy, its current limitations, and future developments.

General organization 220
The SuperflexPy framework has a hierarchical organization with four nested levels: "element", "unit", "node", and "network", collectively referred as "components". These components are shown in Figure 1 and described below: 1. Element (Figure 1a). This level represents the basic model building blocks and is used to create reservoirs, lag functions, and connections. An element can be used to represent an entire 225 https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License. catchment, or, more commonly, a specific process within the catchment.
The reservoir element is described mathematically by ordinary differential equations (ODEs): where S are the state variables (e.g., water storages, substance concentrations, etc.), X are the 230 inputs (e.g., precipitation), Y are the outputs (e.g., streamflow), and S g and Y g are specified constitutive functions (e.g., storage-discharge relationships). In most conceptual models, reservoir elements have a single state variable (representing water storage) however multiple state variables can be accommodated when necessary (e.g., to represent transport).
The lag function element is described mathematically by a convolution integral: 235 where * denotes the convolution operator, X is the input (e.g., water flux), H g is the impulse response function, and T is the time of influence of H g (i.e. the maximum lag). Lag functions are used to represent delays due, for example, to routing.
The connection element joins or splits fluxes from other elements. It has parameters but no 240 states: where C g describes the connectivity between input fluxes and output fluxes, and θ represents connectivity parameters (if any). (Figure 1b). A unit is a collection of multiple connected elements, and is generally intended 245 to implement a lumped catchment model. Multiple reservoir and lag function elements within a unit can be connected to each other, either directly (one-to-one connections), or using connection elements such as splitters and junctions (when a single element is connected to multiple elements). Elements can be combined to build a unit, with the only restriction being that feedback loops between the elements are not allowed. In technical terms, the overall model 250 structure must be an acyclic directional graph. This design restriction is motivated by computational efficiency reasons, as it enables the numerical solution of elements from upstream to downstream. Note that the restriction is not absolute, because it does not preclude feedback between the states within a given element. Hence, if feedbacks are deemed necessary, they can https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.

Unit
handled within an individual element, for example by creating a reservoir element with multiple 255 storages.
3. Node (Figure 1c). A node is a collection of multiple units that operate in parallel. In the context of distributed models, the node can be used to represent a single catchment and the units can be used to represent multiple landscape elements (areas) within the catchment. Each unit within a node is characterized by a weight (e.g. representing its area fraction) specified by the modeler. 260 The weights are used to combine the output fluxes from the units into the total output flux of the node. (Figure 1d). A network connects multiple nodes into a tree structure, and is typically intended to develop a distributed model that generates predictions at internal sub-catchment locations (e.g. to reflect a nested catchment setup). The network routs the fluxes from upstream 265 nodes (leaves of the tree) to the final downstream node (root of the tree). The routing in the river network can be simulated adding delays (lag) to the nodes outputs.

Network
This hierarchical organization makes the effort required to configure SuperflexPy to a new problem proportional to the problem complexity. In particular: • Level 1 is sufficient to create single-element models, e.g., a single-reservoir model or a unit 270 hydrograph model (e.g. Kirchner, 2009); • Level 2 is sufficient to create a lumped model structure, such as GR4J (Perrin et al., 2003) or Hymod (Boyle, 2001); • Level 3 is sufficient create a distributed model that represents spatial heterogeneity but generates predictions only at the catchment outlet (e.g. Gao et al., 2014;Nijzink et al., 2016); 275 • Level 4 is needed only in models that generate predictions at interior points (e.g. Fenicia et al., 2016;Dal Molin et al., 2020b).
Examples of SuperflexPy models implemented at Levels 2 and 4 are given later in Sect. 3.
From a software design prospective, SuperflexPy embraces the object-oriented paradigm (e.g., Meyer, 1988). All framework components are represented by objects that can operate either alone or together, 280 interacting with each other and with external libraries (e.g. for calibration) through defined interfaces.
More details are provided in Sect. 4.2.

A simple illustration of SuperflexPy: creating a new model from existing components
This section illustrates the key steps needed to configure and run a hydrological model using the SuperflexPy framework. The illustration presents a distributed model intended to represent a catchment 285 with 2 HRUs and 3 sub-catchments. The model structure is shown in Figure 1d. Within SuperflexPy, the https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License. entire catchment is represented using a network, the sub-catchments are represented using nodes, and the HRUs are represented using units. The corresponding SuperflexPy code is shown in Figure 2.
For this example, an implementation of the necessary elements with SuperflexPy already exists, therefore the elements only need to be imported. The case where the model structure requires elements for which 290 an implementation is not yet available is considered in Sect. 2.3. Even more complex setups are described in Sect. 3 and in the online documentation (see code availability section).
We start by importing the model components required by the model structure, namely the elements (LinearReservoir and HalfTriangularLag), unit, node, and network. The numerical solvers PegasusPython and ImplicitEulerPython needed to solve the reservoir elements are also 295 imported (more on this in Sect. 4.3). The import operation is shown in lines 1-7.
The imported components are then initialized, which entails specifying the model architecture (connectivity between model components) and the initial values of parameters and states. The initialisation sequence starts from the numerical routines (lines 10-11) and then proceeds from the lowestlevel components (elements) to the highest-level component (network). 300 In detail: L1. Elements are initialized by specifying parameters, states, identifier (id) and, when needed, the numerical solver (lines 14-16).
L2. Units are initialized by specifying the elements that compose them and the identifier (lines 19-20).
The connectivity between elements is defined by conceptualizing the unit as a succession of layers 305 that contain the elements. Further examples on this are given in Sect. 3.
The parameters and states of these elements can be changed after initialization using the methods The model can now be run by calling the get_output method of the highest-level component, as shown 320 on line 42.

Creating new model components with SuperflexPy
We now consider the case where the intended model structure has components beyond those already available in SuperflexPy. New model components can be created by extending existing SuperflexPy components. We anticipate that the SuperflexPy components most likely to require extension are the 325 elements, where new reservoir constitutive functions may be required for new applications. In contrast, it is less likely that unit, node and network functionalities would require extension.
The extension of existing SuperflexPy elements to create new elements relies on the object-oriented paradigm underlying the SuperflexPy software design. The inheritance principle, one of the core concepts of the object-oriented paradigm, allows the user to construct new components by "inheriting" most of the 330 functionalities (methods) from existing classes. Separate implementation is then required only for methods where model differences are to be introduced. This approach reduces substantially the amount of coding required to introduce a new model component. To this end, SuperflexPy provides a library of built-in high-level components that can be easily extended to achieve the desired functionality.
A detailed example of making use of this design is given in Sect. 3.2, which shows how to implement a 335 reservoir with a new storage-discharge relationship.

Examples of building hydrological models using SuperflexPy
This section provides more examples of using SuperflexPy to implement hydrological models, including the use of built-in elements and the creation of new elements. We follow a progression from simple to complex. Section 3.1 shows the implementation of model M4, a lumped model built solely from reservoir 340 elements and used in the original SUPERFLEX case study . Section 3.

Implementing SUPERFLEX configuration M4
M4 is a simple lumped model presented in Kavetski and Fenicia (2011). As shown in Figure 3, M4 comprises two reservoirs connected in series: an "unsaturated" reservoir (UR) intended to represent the partitioning of precipitation between evaporation and runoff, and a "fast" reservoir (FR) intended to represent subsequent streamflow generation mechanisms. , and a portion (UR) Q that is directed to the downstream FR reservoir: where: FR is a power-law reservoir, 360 with the storage-discharge relationship given by ( ) where (FR) k and (FR) α are model parameters.
The inflow (FR) P is given by the outflow from UR, i.e., M4 is a lumped model with multiple elements, and hence can be implemented using SuperflexPy levels L1 and L2 (element and unit, see Section 2.1). Figure 4 shows the code needed to implement M4. Similar to the model described in Sect. 2.2, the two model elements (UR and FR) are already implemented. Hence, the user only needs to import (lines 1-3) and initialize (lines 7-13) the elements together with the numerical https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.
routines. Next, the unit that comprises the two reservoirs is imported (line 4) and initialized (line 15). The 370 input data, namely precipitation and PET time series, are set on line 20. Input data is provided using Numpy arrays. The reading of input data (from text file(s), databases, etc.) is done separately from SuperflexPy, using any suitable Python library or function. In this case, we use Numpy to read from a text file, as shown in lines 17-18.
The model configuration is then complete -line 23 runs the model with given input data to produce the 375 model outputs. The outputs contain streamflow time series in the form of Numpy arrays.

Changing the equations of the fast reservoir in M4
Suppose the modeler wishes to modify model M4 by changing the storage-discharge equation of the fast reservoir given in equation (10) where (FR) k , (FR) α , and (FR) b are model parameters.
An element with this storage-discharge relationship has not been implemented in SuperflexPy yet (as of version 1.2.0). The following sections give two approaches for creating such an element.

Standard approach for creating a new reservoir
The standard approach for creating a new reservoir in SuperflexPy is to define a new class that inherits 385 most of its functionality (methods) from the class ODEsElement. This operation is illustrated in the code snippet in Figure 5. The new class must override the following methods: • __init__: constructor of the class. Its main purpose is to call the constructor of the parent class (lines 5-6) and to point at the method used to calculate the fluxes, here,_fluxes_function_python (see also Sect. 4.3, which shows the benefits of using 390 Numba-optimized methods for calculating the fluxes); • set_input: takes the input fluxes in a predefined order (here, just precipitation) and assigns them a key (line 13) that is then used when setting up and solving the model equations; ODEsElement. Otherwise, in addition to the methods listed above, we would have needed to implement 405 many other methods, e.g., for interfacing with numerical solvers, for setting parameters and states, etc.

Simplified method for creating a new reservoir
The same new reservoir element can be implemented in a simpler way by noting that Note that this simplified implementation is a consequence of the required modification being relatively 415 minor, i.e., a change solely in the constitutive function equation. More complex modifications, such as the inclusion/exclusion of input/output fluxes (e.g. inclusion of evapotranspiration into the PowerReservoir), would require the standard implementation approach described in Sect. 3.2.1.

Implementing a distributed model
This section illustrates the use of SuperflexPy to implement a distributed hydrological model. The 420 example follows the procedure illustrated in Sect. 2.2, creating the more realistic model M02, developed in Dal Molin et al. (2020b) to provide streamflow predictions at 10 sub-catchments of the Thur catchment in Switzerland, see Figure 7a. Each sub-catchment receives its own forcing (precipitation, potential evapotranspiration, and temperature). Two HRU types are defined based on geology: consolidated and unconsolidated formations (Figure 7b). Both HRU types are characterized by the same model structure, 425 which is shown in Figure 8 and represents an enhancement of the structure of model M4 with the following additions: (1) a "snow" reservoir, WR, which controls the partition of incoming precipitation between rainfall and snowfall based on temperature, (2) a lag function between UR and FR, and (3) a https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.
"slow" reservoir, SR, which acts in parallel to FR and is controlled by the same equations as FR but with different parameter values. 430 This example represents a higher degree of complexity compared to the previous examples, both in terms of lumped model structure used for the HRU types (unit) and in terms of introducing a spatial discretization.
Given the spatial organization of the model, nodes are used to represent sub-catchments and units are used to implement HRU types. Note that the sub-catchments may share (one or more) HRU types, which 435 in SuperflexPy translates into the nodes sharing (one or more) units. The network level is used to connect multiple nodes, and enables predictions at internal catchment locations. Figure 10 shows the SuperflexPy representation of the spatial organization shown in Figure 7.
We start by implementing the units. As can be seen in Figure 8, the model structure used to represent the HRUs has elements operating in parallel and, therefore, requires the use of connections. Figure 9 shows 440 how the model structure is "translated" into the SuperflexPy framework. Recall, from Section 2.1, that the connection of elements within a unit must correspond to an acyclic directional graph, i.e., it should not contain feedback loops. Furthermore, elements can be connected only if they belong to two consecutive layers, which implies that "gaps" in the structure must be filled using a transparent element (which outputs the same fluxes it receives as inputs). 445 Comparing Figure 8 with Figure 9, we see how the HRUs structure has been implemented within SuperflexPy. The following implementation aspects are noteworthy: 1. The incoming precipitation is partitioned into rainfall and snowfall. This partitioning is done internally in WR. The SuperflexPy implementation of WR, in fact, takes care of two processes: (i) partitioning of precipitation into rainfall and snowfall; and (ii) simulation of snow processes 450 (accumulation and melting). The output of WR is, logically, the sum of rainfall and snowmelt.
Alternatively, a (new) splitter element could be defined to partition the fluxes between UR (rainfall) and WR (snowfall) based on temperature.
2. WR, as currently implemented, does not receive as input the potential evapotranspiration (PET), which is needed by the downstream element UR. The transfer of the PET to the UR, therefore, is 455 done using the system "upper splitter-upper transparent-upper junction" (Figure 9) that allows to bypass the WR.
3. The parallel part of the structure is composed by two elements on one branch (lag and FR) and only one element on the other branch (SR). To satisfy the requirement of not having "gaps" in the unit structure, a transparent element (lower transparent) is added after. 460 https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.
The code to setup this model is listed in Figure 11. As shown in the simplified example in Sect. 2.2, the user initializes and connects all model components, proceeding sequentially from the lowest level (elements) to the highest level (network). The procedure can be summarized as follows: 1. Lines 10-29: Initialize the elements needed for the lumped model structures used in the HRUs; 2. Lines 32-39: Initialize the units used to represent the HRUs, linking all the elements; 465 3. Lines 42-51: Initialize the nodes used to represent the sub-catchments. Both units are assigned to 9 nodes; the Mosnang sub-catchment contains a single HRU and hence only a single unit is assigned to the corresponding node (line 49).
4. Lines 54-60: Connect the nodes using a network; the topological structure of the network is defined by labeling each node with a unique identifier and specifying the downstream node. 470 The network runs the nodes from upstream to downstream, collects their outputs, and routes them to the outlet. The output of the network is a Python dictionary, with keys given by the node identifiers and values given by the list of Numpy arrays representing the output fluxes.

Implementation details of SuperflexPy
This section presents additional technical details of SuperflexPy. A more detailed description is provided 475 in the model documentation.

Parameters and states
All SuperflexPy components can have parameters and states. Parameters specify component characteristics, whereas states keep track of the component history. States and parameters are set as part of initializing the model components, and can be manipulated using get and set methods provided by 480 the framework at all levels of its hierarchy (see the example in Sect. 2.2).
The parameters can be either constant or variable in time. Constant parameters represent the most common application of hydrological models. Time-variant parameters have been proposed in research applications to represent "deterministic" system variability (e.g. seasonality, Westra et al., 2014) and/or "stochastic" system variability (e.g. Reichert and Mieleitner, 2009). 485

Modular design following the Object-Oriented paradigm
As noted in Sect. 2.3, SuperflexPy embraces the object-oriented paradigm (e.g. Meyer, 1988), which is widely used in general software and is increasingly adopted in scientific software.
• The inheritance principle enables the creation of new classes by extending existing ones. 490 Inheritance reduces drastically the amount of new code that needs to be generated to implement a new model component (an example was provided in Sect. 3.2); • Changes to a class (e.g. a component) and the creation of new classes can be carried out in isolation from the rest of the code, as long as the interfaces between classes are respected; • When creating a model, only the necessary objects need to be initialized and used. This principle 495 makes the model configuration effort roughly proportional to required model complexity, i.e., simple model structures can be constructed from the minimal set of required components. This capability avoids the overhead imposed in frameworks where simple model structures are explicitly constructed as special cases of more complex model structures. The simpler implementation may also reduce computational costs; 500 • Objects retain their history (states), which can be accessed post-run to undertake model analysis and/or subsequent computation; • The modular nature of objects facilitates the development and testing of new code.
These benefits make it easier to achieve clean and maintainable code, which is essential for any practical modelling framework. 505

Numerical solution of ODEs
Reservoir elements are described using ordinary differential equations (ODEs), which are typically solved using numerical time-stepping approximations. There are many such approximations, e.g. Euler methods, Runge-Kutta methods, etc.
SuperflexPy separates the formulation of model equations from the solution of these equations. More 510 specifically, flux equations are defined internally in elements (as shown in the example in Sect. 3.2), while the numerical method is specified externally (to the element) by defining a so-called "numerical approximator". This separation enables the modeller to select the numerical method without making any changes to the model equations.
SuperflexPy provides two built-in numerical approximators, namely the fixed-step implicit and explicit 515 Euler methods. The user can implement additional solvers, either by coding them directly or by interfacing with external ODEs solvers (e.g. from SciPy). As detailed in Sect. 4.4, the choice of numerical implementation, and its compatibility with optimizing compilers, may have a strong impact on the overall computational speed of the model. https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.

Computational efficiency 520
Computational efficiency is a key requirement of a practical modelling framework. Conceptual hydrological models are often used in Monte Carlo uncertainty quantification, which requires 100's -1000's of model runs (even millions in some cases). Model calibration is another common computationally demanding task required by most hydrological models.
The choice of programming languages inevitably requires a trade-off between computational efficiency 525 and ease of use. The choice of Python for SuperflexPy was motivated by the attraction of a flexible and widely used scripting language in conjunction with two efficient numerical libraries: NumPy (Walt et al., 2011) and Numba (Lam et al., 2015). Numpy provides highly efficient arrays for vectorized operations (i.e. elementwise operations between arrays). Numba provides a "just-in-time compiler" that can be used to compile (at runtime) a normal Python method to machine code that interacts efficiently with NumPy 530 arrays.
The combined use of Numpy and Numba is extremely effective when solving ODEs, where the method loops through a vector to perform element-wise operations. The built-in approaches for solving ODEs are compatible with such numerical infrastructure, and therefore enable fast computation times. Note that switching to ODEs solvers that do not take advantage of such libraries might dramatically increase the 535 model runtime.
Numba offers drastic computational speed ups compared to native Python; our experimentation suggests runtime reductions by factors of up to 30. However, a drawback of Numba is the requirement to compile the code at runtime. For a lumped model composed of a few reservoirs, the Numba compilation time is of the order of a few seconds. Therefore, Numba will outperform Python when the simulation is long (e.g. 540 100,000 time steps, corresponding to roughly 12 years of hourly data) and/or when the model needs to be run a large number of times. For example, calibration to observed data (1000's of model runs) takes a few seconds with the Numba implementation compared to a couple of minutes with native Python execution (we here report only the runtime of the model itself, and exclude the runtime of the calibration tool procedures; more details on benchmarking in the documentation). 545

Ability to represent multiple fluxes and states
SuperflexPy can operate with multiple fluxes and state variables. In particular, connection elements, units, nodes, and the network are designed to deal with an arbitrary large number of fluxes.
This generality is intended to support the extension of SuperflexPy to model the transport of chemical substances contained in the water (e.g., Fenicia et al., 2010;Ammann et al., 2020). Representing fluxes of 550 chemical substances requires state variables and fluxes in addition to those corresponding to water https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.
storages and water fluxes. The calculation of such fluxes also requires additional equations and parameters. The available examples in SuperflexPy do not include transport processes. However, the framework architecture foresee this need, and has been designed to be readily extended to accommodate such processes. For example, elements (e.g. a reservoir) can be created to implement the functionality 555 required to simulate such fluxes (e.g. specifying the governing equations for substance transport, etc).

Defining the right complexity for a flexible model implementation
Achieving the envisaged flexibility of the SUPERFLEX framework in practical software is challenging, because flexible modelling functionality may come at the expense of ease of use and computational cost. 560 The challenge in designing SuperflexPy has been to determine an appropriate level of abstraction for typical conceptual model applications. On one hand, high level of mathematical generality and abstraction offers the most flexibility, but risks losing a clear hydrological interpretation and may increase user effort in customizing the model. On the other hand, a framework that is over-restricted in terms of component behaviour may be easy to manage (as the number of modelling options is low), but it may not fulfil the 565 promise of flexible models and may result in a limited range of application (e.g. limiting to lumped configurations only).
Conceptual models vary significantly in terms of complexity (e.g. from single bucket models to distributed models, from modelling streamflow alone to modelling water isotopes and/or other chemicals, etc.). In order to accommodate this range of potential model complexity in a flexible and practical way, 570 we have organized the SuperflexPy software into four hierarchical levels, namely element, unit, node, network. As shown in Sect. 2.1 and in the examples of Sect. 3, these levels map to many types of conceptual modelling applications, from a single element (e.g. a reservoir), to a lumped model (typically composed of several elements, such as a combination of reservoirs), to a composition of lumped models, designed to provide prediction at a single outlet (e.g. a catchment with several HRUs, characterized by 575 lumped models), and eventually to a distributed model capable of making predictions at multiple internal sub-catchments. An important outcome of this design choice is that the model configuration effort by the user becomes proportional to the required model complexity, so that simple models are much easier to configure.
The flexibility in customizing SuperflexPy elements is enhanced through its Object-Oriented design. As 580 shown in Section 3.2, new components can be built by inheriting most methods from existing components and specifying only the required new features of interest. The potential downsides of using a scripted language Python in terms of computational speed are mitigated by the ability to use the Numba package.

Current limitations in model structure specification
As part of balancing the flexibility, ease of use, and computational performance of SuperflexPy, some 585 restrictions have been imposed on the connectivity between model components.
The first restriction is that the elements within the unit must represent an acyclic directional graph, with no feedback loops from downstream to upstream elements (Sect. 2.1). This restriction enables the numerical solvers to proceed in a single pass from upstream to downstream and improves the computational performance of the framework. It also simplifies the specification of its structure. The 590 restriction on internal model feedbacks is not expected to be overly limiting when developing conceptual hydrological models, as the fluxes in these models typically flow only in one direction (e.g. in the model M4 the water flows only from UR to FR and not vice versa). An example where internal model feedbacks may be required is given by the bidirectional interaction between surface water and groundwater in the hyporheic zone, where the direction of the fluxes changes depending on the state of the two components. 595 Such interactions can still be modelled in SuperflexPy by introducing elements that embed such feedbacks internally. For example, the hyporheic zone, can be represented using a two-state reservoir with interacting states (e.g., Seibert et al., 2003). In other words, the restriction on model feedbacks applies to interactions between elements but not to the internal structure within an element.
The second restriction regards the topology of a network, which must represent a tree where any given 600 node can connect and transfer fluxes to only a single downstream node (Sect. 2.1). This requirement is met by typical conceptual distributed models, as discussed in the introduction and illustrated in Sect. 3.3.
However, fully integrated distributed models, such as Parflow or Mike She, do consider mutual dependencies between spatial elements, leading for example to 2D or 3D groundwater flows. Such configurations are considered beyond the scope of conceptual distributed models, and therefore are 605 currently not supported in SuperflexPy.

Current usage and future developments
SuperflexPy is easy to install and run; it is written in pure Python and its dependencies are limited to the packages NumPy and Numba (Sect. 4.4). Installation can be done directly using the package installer for Python (pip) without the need of installing external libraries. We stress that SuperflexPy is not a wrapper 610 of earlier SUPERFLEX code but offers a completely new implementation that is not constrained by choices taken in the earlier code versions.
SuperflexPy has already been used for research applications. Jansen et al. (2020) performed a "model mimicry" study where similarities and differences within the HBV "family" models were investigated.
SuperflexPy was used to construct a compare a set of HBV-like models, assessing, among other things, 615 https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License. the behavior of single components and the impact of numerical implementation. A list of publications using SuperflexPy is maintained on the documentation website.
In terms of future developments, we hope that SuperflexPy offers the hydrological community a new tool for research work and practical applications. Further SuperflexPy developments are likely to follow from such work and collaborations, including: (1) expansion of the library of model components beyond the 620 ones here presented (as shown in the example in Sect. 3.2), and (2) more fundamental changes in response to future model applications. It is important to highlight that SuperflexPy can be used to create and combine new model components, thereby enabling experimentation with new model structures and general conceptualizations. The framework, therefore, is not limited to components and structures taken from existing models. Such collection, however, may be produced by storing the configurations that allow 625 reproducing such models. The model documentation already contains a small sample of such configurations, which may grow as new users share their implementations with the community. In order to facilitate the use of the code, the code is accessible on GitHub with license LGPL-3.0 and distributed using the Python package installer PyPI (refer to the section at the end of the paper on code availability).
The documentation provides a guide on how interested colleagues can contribute to the framework. 630

Conclusions
SuperflexPy is a new Python flexible modelling framework for building conceptual models ranging from lumped to distributed configurations. The framework offers detailed control over each aspect of model configuration, and caters to a wide range of typical conceptual model applications. In order to facilitate the model building process, the framework is organized using four hierarchical levels, namely element, 635 unit, node and network, which correspond to conceptual model setups of increasing levels of complexity, namely a single element (e.g. a reservoir), a lumped model (e.g. a collection of interconnected reservoirs), a collections of lumped models designed to provide prediction at a single outlet, and a distributed model designed to provide predictions at internal sub-catchments. As the construction of a model that requires a certain hierarchical level does not require specifying the levels above it, the model configuration effort is 640 proportional to the complexity of the application. The framework supports multiple states and fluxes in each component, and foresees an extension to water quality applications.
SuperflexPy builds on the experience of the authors and their colleagues in a series of hydrological case studies using the SUPERFLEX principles. SuperflexPy offers an open source implementation of the SUPERFLEX principles. In order to facilitate usability and further developments, we focused on several 645 aspects of the model code, including ease of use and interfacing, availability, amenability of extensions, and computational efficiency. https://doi.org/10.5194/gmd-2020-409 Preprint. Discussion started: 16 December 2020 c Author(s) 2020. CC BY 4.0 License.
We presented two examples illustrating common use cases of the SuperflexPy framework, including the implementation of a simple lumped reservoir model, and the implementation of a distributed model to represent a system of multiple sub-catchments and HRUs. It is hoped that the framework will contribute 650 to ongoing efforts in the hydrological modelling community to develop more robust and representative models. The framework is open source, available with license LGPL-3.0 on GitHub.

Code availability
The organization of SuperflexPy as a software project is shown in Figure 12 using the blue and yellow boxes. The nodes, used to represent the sub-catchments, are shown using the green dashed boxes. The group of nodes connected together (green arrows) creates a network.