MODEL CHECKING FOR REAL-TIME ATTACK DETECTION IN WATER DISTRIBUTION SYSTEMS

. Water distribution systems represent critical infrastructures. These architectures are really critical, and irregular behaviour can be reflected in human safety. As a matter of fact, an attacker obtaining control of such an architecture is able to perpetrate a plethora of damages, both to the infrastructure and people. In this paper, we propose an approach to identify irregular behaviours focused on water distribution systems. The designed approach considers a formal verification environment. The logs retrieved from water distribution systems are parsed into a formal model and, by exploiting timed temporal logic, we characterize the behaviour of a water distribution system while an attack is happening. The evaluation, referred to a water distribution system, confirmed the effectiveness of the designed approach in the identification of three different irregular behaviours.

for Collaborative Enterprises (WETICE). With respect to the paper in [19] below, we depict the contributions of this paper: -we propose a property to identify the reply irregular behaviour; -we extend the experiment to validate better the underflow and the overflow properties presented in [19]; -we introduce a new property for the identification of a new anomaly, i.e., the reply irregular behaviour.
The work proceeds as follows: the next section presents background concepts about model checking timed automata, Section 3 describes the designed approach, Section 4 discusses the performed experiment aimed to evaluate the approach and, finally, in Section 5 conclusions and future research directions are drawn.
2. Background. In this section, fundamental concepts related to the formal methods technique adopted by the proposed method are provided, i.e., the model checking timed automata.
The toolchain to apply the model-checking technique is composed of a formal model and a temporal logic: both of them are described in this section.
We consider timed automata, proposed by Alur and Dill [20,21], as a formal verification technique for real-time systems as, for instance, the water distribution systems. In a nutshell, a timed automaton is composed of a classical finite automaton able to manage clocks, evolving continuously and synchronously with respect to absolute time. Each transition of a timed automaton is labelled by a guard or constraint overclock values, indicating the time in which the transition can be fired, and a set of clocks to be reset when the transition is fired. Each location is constrained by an invariant, which restricts the possible values of the clocks for being in the state, which can then enforce a transition to be taken [22][23][24].
In this paper, we model a sequence of logs gathered from a water distribution systems as a network of timed automata, i.e., a finite-state machine extended with clock variables [21]. The model is extended with bounded discrete variables that are part of the state. The state of the system is defined by the location of all automata, the clock values and the values of the discrete variables. Automaton may fire an edge (i.e., perform a transition) separately or synchronise with another automaton with the aim to lead a new state.
Below we provide the definition of the syntax and semantics for the basic timed automata. We exploit the following notation: C is a set of clocks and B(C) is the set of conjunctions over simple conditions of the form x ▷◁ c or x − y ▷◁ c, where x, y ∈ C, c ∈ N and ▷◁ ∈ { <, ⩽, =, >, ⩾ }. A timed automaton is a finite directed graph annotated with conditions over and resets of non-negative real-valued clocks [21]. Definition 1. Timed Automaton. A timed automaton is a tuple (L, l 0 , C, A, E, I) where L is a set of locations, l 0 ∈ L is the initial location, C is the set of clocks, A is a set of actions, co-actions and internal τ -action, E ⊆ L × A × B(C) × 2 c × L is a set of edges between locations with an action, a guard and a set of clocks to be reset, and I : L → B(C) assigns invariants to locations.
In the following, we define the semantics of a timed automaton. A clock valuation is a function u : C → R ⩾0 from the set of clocks to the non-negative reals. Let R C be the set of all clock valuation. Let u 0 (x) = 0 for all x ∈ C. We will abuse the notation by considering guards and invariants as sets of clock valuations, writing u ∈ I(l) to mean that u satisfies I(l).
Definition 2. Semantics of Timed Automaton. Let = (L, l 0 , C, A, E, I) be a timed automaton. The semantics is defined as a labelled transition system ⟨S, s 0 , →⟩, where S ⊆ L × R C is the set of states, s 0 = (l 0 , u 0 ) is the initial state, and →⊆ S × (R ⩾0 ∪ A) × S is the transition relation such that: where for d ∈ R ⩾0 , u + d maps each clock x in C to the value u(x) + d, and [r → 0]u denotes the clock valuation which maps each clock in r to 0 and agrees with u over C\r. Timed automata are often composed into a network of timed automata over a common set of clocks and actions, consisting of n timed automata i = (L i , l 0 i , C, A, E i , I i ), 1 ⩽ i ⩽ n. A location vector is a vectorl = (l 0 1 , . . . , l 0 n ). We compose the invariant functions into a common function over location vectors I(l) = ∧ i I i (l i ). We writel [l ′ i /l i ] to denote the vector where the i-th element l i ofl is replaced by l ′ i . In the following, we define the semantics of a network of timed automata.
Definition 3. Semantics of a network of Timed Automata: Let = (L, l 0 i , C, A, E i , I i ) be a network of n timed automata. Letl 0 = (l 0 1 , . . . , l 0 n ) be the initial location vector. The semantics is defined as a transition system ⟨S, s 0 , →⟩, where S = (L 1 x . . . xL n ) × R C is the set of states, s 0 = (l 0 , u 0 ) is the initial state, and →⊆ S × S is the transition relation defined by: Modelling languages usually extend timed automata with the following additional features: binary synchronisation: channels are declared as chan c. An edge labelled with c! synchronises with another labelled c?. A synchronisation pair is chosen non-deterministically if several combinations are enabled; broadcast channels: are declared as broadcast chan c. In a broadcast synchronisation one sender c! can synchronise with an arbitrary number of receivers c?. Any receiver that can synchronise in the current state must do so. If there are no receivers, then the sender can still execute the c! action, i.e. broadcast sending is never blocking; initialisers: are used to initialise integer variables and arrays of integer variables. For instance, int i := 2; or int i [3] := {1, 2, 3}.
Once defined the formal model, we need to verify the model with regard to a requirement specification. Similarly to the model, the requirement specification must be expressed in a formally well-defined language. Several such logics exist in the scientific literature. In this paper, we consider the Timed Computational Temporal Logic (TCTL) [21,25], which extends the classical untimed branching-time logic CTL [26] with time constraints on modalities.
The TCTL syntax is defined by the following grammar: where a ∈ AP (we denote with AP a set of atomic propositions), and I is an interval of R + with integral bounds. There are two possible semantics for TCTL, one which is said continuous, and the other one which is more discrete and is said pointwise. We consider the second one, i.e., the pointwise semantic (Table 1).
Where ϱ[π] is the state of ϱ at position π, and duration ϱ ⩽ π is the prefix of ϱ ending at position π, and duration(ϱ ⩽ π ) is the sum of all delays along ϱ up to position π.
In the pointwise semantics, a position in a run: ϱ = s 0 τ 1 , e 1 s 1 τ 2 , e 2 s 2 . . . s n−1 τ n , e n s n , is an integer i and the corresponding state s i . In this semantics, formulas are checked only right after a discrete action has been done. Sometimes, the pointwise semantics is given in terms of actions and timed words, but it does not change anything. As usually in CTL, TCTL: ≡ a ∨ ¬ a standing for true, ≡ ¬ standing for false, the implication ϕ ⇝ ψ ≡ (¬ ϕ ∨ ψ), the eventual operator F I ϕ ≡ tt U I ϕ and the global operator G I ϕ ≡ ¬ (F I ¬ ϕ).
Formulae can be classified into reachability, safety and liveness. Reachability properties: they are looking if whether a given state formula, φ, possibly can be satisfied by any reachable state. Reachability properties are often used while designing a model to perform sanity checks. We express that some state satisfying φ should be reachable using the path formula E φ.
Safety properties are expressed in the following form: "something bad will never happen". These properties are usually formulated positively, for example, something good is invariantly true. For instance, let φ a state formula, we express that φ should be true in all states that are reachable with the path formulae A[ ] φ.
Liveness properties are of the following form: something will eventually happen. Liveness is expressed with the path formula A φ meaning φ is eventually satisfied. A useful form is the leads to or response property, written ψ ⇝ φ which is read as whenever ψ is satisfied, then eventually φ will be satisfied.   Once the model and temporal logic properties are defined, we need something enabling us to check whether the timed automata network (i.e., the formal model) satisfies the defined properties. To this aim, we consider formal verification, a system process exploiting mathematical reasoning to verify whether a system under analysis (i.e., the model) satisfies some requirements (i.e., the timed temporal properties).
Several verification techniques have been proposed in the last years. In this paper, we resort to model checking [24,27].
In the model checking technique the properties are formulated in temporal logic: each property is evaluated against the system. The model checker accepts as input a model and a property, it returns "true" whether the system satisfies the formula and "false" otherwise. The performed check is an exhaustive state space search that is guaranteed to terminate since the model is finite. In this paper, we consider as model checker UPPAAL 1 [22,23,28], an integrated tool environment for modeling, validation and verification of real-time systems modeled as timed automata networks. Thus, the syntax of UPPAAL expression is given in Table 2.
Expressions are used with the following labels: guard: a particular expression satisfying the following conditions: (i) it is side-effect free; (ii) it evaluates to a boolean; (iii) only clocks, integer variables, and constants are referenced (or arrays of these types); (iv) clocks and clock differences are only compared to integer expressions; (v) guards overclocks are essentially conjunctions (disjunctions are allowed over integer conditions); synchronisation: a synchronisation label is either on the form Expression! or Expression? or is an empty label. The expression must be side-effect free, evaluate to a channel, and only refer to integers, constants and channels; assignment: an assignment label is a comma-separated list of expressions with a side-effect; expressions must only refer to clocks, integer variables, and constants and only assign integer values to clocks. invariant: an expression that satisfies the following conditions: it is a side-effect free; only clock, integer variables, and constants are referenced; it is a conjunction of conditions of the form x < e or x ⩽ e where x is a clock reference and e evaluates to an integer.
3. The Designed approach. In the designed approach, we consider the cistern level gauging as features. In details, if the water distribution system is formed by two cisterns, i.e., two features describing the water gauging in the two cisterns are considered.
The designed approach consists of two main steps: the Formal Model Creation ( Figure 5) and the Formal Model Verification (Figure 8).
From the water distribution system under analysis and the technical report, a technician marks the specific day log as irregular. The features gathered from the cistern levels and the label to mark a specific trace as irregular are the input for the one day log stored in CSV files, containing the cistern gauging for one day at a fixed range time (1 hour).
The Discretisation is aimed to discretise each cistern level feature. The numeric values are split into 3 different ranges [29]. Furthermore, each discretised feature previously gathered is converted into a timed model. The discretisation process is aimed to convert continuous values into discrete ones. A plethora of methods were proposed by the research community for numeric values discretisation: in our approach, we exploit one discussed by researchers in [29]. The idea behind this method divides the numeric features into three intervals by exploiting the equal-width partitioning that basically divides the values of a certain numeric value into three equal-size intervals. In our case, the cistern level continuous values are divided into one of three different classes (Up, Basal, and Low). The partitioning is applied to all features (i.e., the cisterns). Moreover, once obtained the discrete values from the discretisation process, from each feature is generated a timed automaton (i.e., formal model in Figure 5).
With the aim to better explain how the timed automaton is built, let us consider the following example. In particular, in Table 3 we represent a fragment of a discretised feature.
With the first column (i.e., Time in Table 3) we are referring to the interval time (in the fragment 1 ⩽ t ⩽ 6), while with the F 1 and F 2 columns, we are referring to the two considered features (i.e., the cistern levels). By considering the fragment in Table 3, in the t 3 time interval, the F 1 reaches the U p value, while the F 2 feature reaches the Low value. Once obtained the discretised values for the features (i.e., the cistern values), it is possible to generate the formal model. In detail, we build a timed automata network: in a nutshell, for each discretised feature, we generated an automaton.
Basal Up t 5 Low Basal t 6 Up Basal By considering the discretised fragment, we shown in Table 3 the timed automata network that is built in the following way: if the same discrete value is repeated in a consecutive time interval, the automaton related to the feature under analysis exhibits a loop: for instance, the automaton built from the F 1 discretised feature contains a loop related to the t 1 , t 2 and t 3 time intervals (because the Up discretised values are repeated three times). Moreover, the resulting automaton for the F 2 feature shows two loops: the first loop is related to the t 2 and t 3 time intervals (Low is the repeated discrete value), while the second loop is related to the t 5 and t 6 time intervals, and in this last case the repeated discrete value is Basal).
In order to exit from the loops, a guard is exploited. Differently, with regard to the entering condition, an invariant is considered. Moreover, for each automaton, we consider a number of clocks equal to two: the first clock (x) to ensure the entering condition into the loop and the second one (y) to ensure the exit condition. Moreover, each automaton is responsible for storing the count related to the respective Up, Basal and Low values.
The values for the F x automaton are labelled with a subscript x ∈ {1, 2} (we consider two features). The s channel permits the automata synchronisation and, for this reason, is not stored locally. The aim of the channels is to ensure the automata network progression. This mechanism basically is a hand-shaking synchronization: two processes take a transition at the same time, the first one will have an s!, while the other an s?, in order to ensure the synchronization. As a matter of fact, one sender event s! is able to synchronise with a number of s? receiver events. We resort to channels in order to avert incoherence from the discretised feature values and the interval times. For instance, the automaton related to the first feature is indicated with F 1 and its variables are u 1 , b 1 and l 1 respectively for Up, Basal and Low values. As a consequence, the automaton for the second variable is F 1 , and its variables are u 2 , b 2 and l 2 , respectively for Up, Basal and Low values of this second automaton. Figures 6 and 7 indicate the model, respectively, gathered from the F 1 and F 2 discretisation. In both the figures, we indicate with u the Up value, with b the Basal value, and with l the Low value, all the three values are present in Table 3.
Below we provide a brief explanation about the models represented in Figures 6 and 7: the F 1 automaton is iterating in the loop (node 1 in Figure 6) for three-time intervals (i.e., y 1 < 3), while the F 2 automaton after the increment of the u 1 variable (node 1 in Figure 7) is iterating for two-time intervals (i.e., y 2 < 2). Subsequently, the F 1 automaton does not exhibit any loops, while the F 2 automaton is iterating for two-time intervals in the loop in node 3 in Figure  7, and then it continues with the last node. Once generated the formal model, the Formal Model Verification step (Figure 8) is aimed at checking if the modeled water distribution system is violated. The Formal Model Verification receives as input a formal model (formal model) built in the previous step and a set of properties. The set of logic properties is checked against the formal model generated (formal verification environment) using the UPPAAL (an acronym based on a combination of UPPsala and AALborg universities) formal verification environment 2 . UP-PAAL represents an integrated tool environment for modeling, validation and verification of real-time systems modeled as networks of timed automata.
4. Experimental evaluation. From the physical layout point of view, the following generic scenario is considered: the analysed SCADA water distribution system is composed by a generic water distribution system operator where, recently, has introduced a new technology to enable remote data collection from sensors and remote control of actuators. After the technology was introduced, anomalous levels in two tanks were observed; for instance, water overflow in one tank occurred. By looking for the causes, experts domain suspect potential cyberattacks. In particular, they consider some kind of malicious behaviours [30,31] aimed to activate and deactivate the actuators.
The dataset 3 exploited in the experimentation of the designed approach contains one-day logs from a water distribution system composed from the water levels of the two cisterns (i.e., ct1 and ct2).
We have hourly gauging of the water under regular operating conditions, and when an overflow (OF ), an underflow (U F ) or a reply (R) irregular behaviours targeting, respectively, the ct1 and ct2 cisterns are happening. Logs referred to 30 days of gauging are considered: ten marked with the OF irregular behaviour on ct1, ten with the U F irregular behaviour on ct2 and the last ten days with reply (R) irregular behaviour on both the ct1 and ct2 cisterns (irregular behaviours have happened each day), while the remaining ten days without irregular behaviours [5,32,33]. For each day log a formal model is built. In total, we consider 40 days.
To summarise, in the experiment, we consider 40 days: 30 days related to irregular behaviours and 10 days with legitimate behaviours. The 30 days contain the following irregular behaviours: 10 days with overflow attack each day, 10 days with underflow attack each day and 10 days with reply attack each day.
4.1. The Properties. Table 4 indicates the OF , U F and R irregular behaviour identification properties.
The φ property, referred to the OF irregular behaviour, is able to verify if an OF is in progress (i.e., when the ct1 level exhibits an up value at least 11 times). The U F irregular behaviour, described by the χ property, is able to verify if an U F is in progress (i.e., when the ct2 level exhibits a low value at least 9 times). The last property, i.e., ψ (for the R irregular behaviour identification), is aimed to verify if both the ct1 and the ct2 cisterns present the same value for more than 12 times. We want to verify if the φ, the χ and the ψ properties, possibly, can be satisfied by a reachable state. We express that some state satisfying φ should be reachable using the path property Eφ. Similar considerations can be done for the χ and ψ properties. Furthermore, provide the property expressing that can happen an OF , or U F or a R irregular behaviours (i.e., Eφ ∨ χ ∨ ψ).
The properties were formulated with the help of domain experts. As a matter of fact, the properties are aimed to explain the domain experts knowledge. In particular, the domain experts suggest that whether ct1 exhibits for a number of times equal or greater than 11 an increasing value, this is reflecting of an overflow attack in progress, while if ct2 exhibits a value that is decreasing for a number of times equal or greater than 9, this is reflecting in an underflow attack in progress. With regard to the reply attack, expert domains suggested that whether both ct1 and ct2 are showing basal values for a number of times equal or greater than 12 this is symptomatic for this kind of attack in progress.

The Experiment.
The property verification outcomes are indicated in Table 5.
With L x we mark the log of the water cisterns for a single log, where L ∈ {irregular, regular}, x identifies the log under analysis. The dataset comprises 40-day log, 10 exhibiting normal behaviour, while the remaining 30 are afflicted by OF , U F and R irregular behaviours. The column (φ) identifies the models verified as true (✓ in Table 5) or false (✗ in Table 5) to the OF property, while the column (χ) identifies the models verified as true or false to the U F property.
From the formal verification environment output indicated in Table 5 we can state the φ and the χ are able to identify all the models afflicted by the OF irregular behaviour (i.e., irregular 1 o and irregular 10 o ). The χ property is able to identify all the models afflicted by the U F irregular behaviour (i.e., irregular 1 u , irregular 2 u , irregular 3 u , irregular 4 u , irregular 5 u , irregular 6 u , irregular 7 u , irregular 8 u , irregular 9 u and irregular 10 u ). The ψ property is able to identify all the models afflicted by the R irregular behaviour (i.e., irregular 1 r , irregular 2 r , irregular 3 r , irregular 4 r , irregular 5 r , irregular 6 r , irregular 7 r , irregular 8 r , irregular 9 r and irregular 10 r ). Furthermore, all the models without irregular behaviour (i.e., regular 1 , regular 2 , regular 3 , regular 4 , regular 5 , regular 6 , regular 7 , regular 8 , regular 9 and regular 10 ) are marked as false by φ, χ and ψ properties.
Moreover, we consider four different metrics to evaluate the performance of the proposed approach: Precision, Recall, F-Measure and Accuracy.
The precision has been computed as the proportion of the observations that truly belong to investigated logs among all those which were assigned to the specific attack. It is the ratio of the number of relevant records retrieved to the total number of irrelevant and relevant records retrieved: where tp indicates the number of true positives (for instance, whether w e are evaluating the φ formula, this value represents the number of one-day logs whose related overflow attack model is correctly labelled as true by the formal verification environment) and fp indicates the number of false positives (for instance, whether we are evaluating the φ formula, this value represents the number of models whose related legitimate model is wrongly labelled as true by the formal verification environment).
The recall has been computed as the proportion of attacks that were assigned to a given class among all the attacks that truly belong to the class. It is the ratio of the number of relevant records retrieved to the total number of relevant records: where tp indicates the number of true positives and fn indicates the number of false negatives (for instance, whether we are evaluating the φ formula, this value represents the number of one-day log whose related overflow attack model is wrongly labelled as false by the formal verification environment). The F-Measure is a measure of a test's accuracy. This score can be interpreted as a weighted average of the precision and recall: where Precision and Recall are obtained by following the formulae above explained. The Accuracy is the fraction of the model correctly identified, and it is computed as the sum of true positives and negatives divided all the evaluated models: lclAccuracy = tp + tn tp + f n + f p + tn , where tn indicates the number of true negatives (for instance, whether we are evaluating the φ formula, this value represents the number of one-day logs whose related legitimate model is correctly labelled as false by the formal verification environment). By the metrics computation, we reach a value equal to 1 for all metrics we have considered (i.e., Precision, Recall, F-Measure and Accuracy). This is symptomatic that the designed approach is able to correctly identify the several irregular behaviours we considered in the experiment with no false positive.

Conclusion and future work.
The risks of an attack on SCADA systems are concrete and real: these are very sensitive objectives given the importance of the processes they govern and the impact that a disservice can cause on the community. Precisely, because of their structure distributed over the territory, hitting a SCADA system can really affect hundreds or thousands of other sites or plants, which in turn become unmanageable and therefore dangerous. There are countless organizations and goals that move interests in carrying out or making increasingly sophisticated actions and attack strategies, ranging from fraud and physical to creating disservices.
For this reason, to mitigate the irregular behaviours in critical infrastructure and, in particular, in SCADA water distribution systems, in this paper, we design a timed automata-based approach to identify OF , U F and R irregular behaviours on the water distribution system. We propose formal model logs obtained from SCADA systems in terms of a timed automata network by exploiting the UPPAAL formal verification environment.
An accuracy equal to 1 is reached, symptomatic of the effectiveness of the proposed approach in attack detection in water distribution systems.
The proposed method can be adopted for real-time irregular behaviour detection in critical systems while the monitored SCADA system is under attack. As a matter of fact, in critical contexts, it is of fundamental importance to identify a threat when the latter is in progress in order to minimize data and put the system in a position to work safely. Moreover, considering that the proposed approach exploits formal methods and temporal logic formula, it is possible to understand the reason why a certain property is resulting false using the counterexample. Although the principal aim of the counterexample is to assist the designer in finding the source of the error in complex systems design, the counterexample can be exploited for other purposes. For instance, in the context of the proposed method can be used to understand how many time intervals the property will become true, thus providing a sort of probability of risk that irregular behavior may occur. This aspect can be of interest because it can be considered to forecast future irregular behaviours, thus taking measures before they can possibly take place.
As future work, we plan to evaluate the proposed of other SCADA critical systems (for instance, relating to oil/gas or power management and distribution) with the aim to validate the proposed method in another critical environment. Moreover, we plan to evaluate the proposed method with systems with much faster dynamics (for instance, power distribution/generation), i.e., with significantly larger models.