VISIR-I: small vessels – least-time nautical routes using wave forecasts

<jats:p>Abstract. A new numerical model for the on-demand computation of optimal ship routes based on sea-state forecasts has been developed. The model, named VISIR (discoVerIng Safe and effIcient Routes) is designed to support decision-makers when planning a marine voyage. The first version of the system, VISIR-I, considers medium and small motor vessels with lengths of up to a few tens of metres and a displacement hull. The model is comprised of three components: a route optimization algorithm, a mechanical model of the ship, and a processor of the environmental fields. The optimization algorithm is based on a graph-search method with time-dependent edge weights. The algorithm is also able to compute a voluntary ship speed reduction. The ship model accounts for calm water and added wave resistance by making use of just the principal particulars of the vessel as input parameters. It also checks the optimal route for parametric roll, pure loss of stability, and surfriding/broaching-to hazard conditions. The processor of the environmental fields employs significant wave height, wave spectrum peak period, and wave direction forecast fields as input. The topological issues of coastal navigation (islands, peninsulas, narrow passages) are addressed. Examples of VISIR-I routes in the Mediterranean Sea are provided. The optimal route may be longer in terms of miles sailed and yet it is faster and safer than the geodetic route between the same departure and arrival locations. Time savings up to 2.7 % and route lengthening up to 3.2 % are found for the case studies analysed. However, there is no upper bound for the magnitude of the changes of such route metrics, which especially in case of extreme sea states can be much greater. Route diversions result from the safety constraints and the fact that the algorithm takes into account the full temporal evolution and spatial variability of the environmental fields.                     </jats:p>

section_content: The operational availability of high spatial and temporal resolution forecasts, for weather, sea state, and oceanographic variables paves the way to a realm of downstream services, which are increasingly closer to end-user needs (Ryder, 2007). Such services may support the decision-making process in critical situations where knowledge of the present and predicted environmental state is key to avoiding casualties or to making savings in terms of time, economic cost, or environmental impact. 

VISIR [vi'zi:r] 1 is a model 2 and an operational system 3 for the on-demand computation of safe and efficient ship routes based on sea-state forecasts. In its present version, VISIR-I, medium and small motor vessels with displacement hulls are considered, such as fishing vessels (e. g. sein-G. Mannarini et al. : VISIR-I: least-time nautical routes ers, trawlers), towboats and fireboats, service boats (crew and supply boats), short trip coastal freighters, displacement hull yachts and pleasure crafts, and small ferry boats. 

The aim of this paper is to lay a sound scientific foundation of VISIR-I, including all its main components: the optimization algorithm, the ship model, and the processor of the environmental fields. 

After reviewing the literature in Sect. 1. 1 and summarizing our original contribution in Sect. 1. 2, the solution devised for VISIR-I is presented in detail in Sect. 2. Examples of optimal routes in the Mediterranean Sea (Sect. 3) precede the conclusions, which are drawn in Sect. 4.
section_name: Introduction

section_content: The main mathematical schemes available in the literature to solve ship routing problems are reviewed in the following. 

Initially devised as a manual tool for navigators, the isochrone method is based on the idea of building an envelope of positions attainable by a vessel at a given time lag after departure. This envelope is called an "isochrone". In the work by Hagiwara (1989), a detailed algorithm is provided, describing how to generate the isochrones and how to use them for constructing a least-time route. Space and course discretization in the vicinity of the rhumb line between departure and arrival locations are performed. At each progress stage, the course leading to the maximum spatial advancement from the origin is considered. When an isochrone gets close enough to the destination, the optimal route is recovered by a backtracking procedure. No proof of the time optimality of the resulting route is provided. Hagiwara's modified isochrone method is the basis for the fuel optimization method proposed by Klompstra et al. (1992). Here, each stage is represented by a two-dimensional position and time. Instead of isochrones or time fronts, energy fronts or "isopones" are computed, being the attainable regions for a given expenditure on fuel. Szlapczynska and Smierzchalski (2007) review several variants of the isochrone method, highlighting their weaknesses, such as limitations in the form of ship speed characteristics and in dealing with landmasses, especially in the vicinity of narrow straits. The authors propose a solution to the latter issue, by screening all route portions intersecting the landmass. 

The variational approach involves searching for trajectories making an objective functional stationary, such as total time of navigation or operational cost, given a set of constraints. The search is achieved by varying the parameters controlling the trajectory. This approach is equivalent to solving the Euler-Lagrange equation. In Hamilton (1962), leasttime ship routes are computed by varying the ship's course, under the assumption that the environmental field is static and thus vessel speed does not explicitly depend on time. 

The time-dependent problem instead can be addressed through the technique of optimal control (Pontriagin et al., 1962). With this method, the dynamic system (the vessel) is controlled by a time-dependent input function (typically engine thrust and rudder angle), allowing the objective function to be minimized. Optimal control is formulated in terms of a set of necessary conditions (Luenberger, 1979). Applications of optimal control to ship routing problems are found in Bijlsma (1975), Perakis and Papadakis (1989) and Techy (2011). Least-time transatlantic routes are computed by Bijlsma (1975). There, significant wave height and wave direction fields from 12-hourly forecasts are assumed to determine vessel speed, while the sole control variable is vessel course. The method can account for prohibited courses due to dynamic reasons (e. g. rolling). However, specific geometrical conditions on the vessel speed characteristics have to hold for the method to work. Furthermore, due to topological issues, there are unreachable regions of the ocean, and the method involves guessing the initial vessel course, which may hinder the implementation in an automated system. The approach by Perakis and Papadakis (1989) accounts for a delayed departure time and for passage through an intermediate location (point-constrained problem). Local optimality conditions ("broken extremals") are found at the boundaries of spatial sub-domains. The optimal ship power setting is found to always take the maximum value possible. The results hold under the assumption that the ship speed characteristics depend on engine throttle as a multiplicative factor. Another limitation of this approach is that the computed extremal trajectory is not guaranteed to lead to a minimum of the objective function. In Techy (2011) the author reports on a vessel moving with constant velocity with respect to water in presence of currents ("Zermelo's problem"). The optimal trajectory is analysed as a function of flow divergence and vorticity, finding the optimal steering policy in a point-symmetric, time-varying flow field. In addition, a geometrical interpretation of Pontriagin's principle is provided. However, to deliver a unique solution, the method requires the hypothesis that the domain of maneuverability of the ship is convex. 

The work by Lolla et al. (2014) is based on the computation of the reachability front of a vehicle with an internal propulsion system, subject to a time-dependent ocean flow. The front is implicitly defined through a level set, and its evolution satisfies a specific solution of a Hamilton-Jacobi equation. The optimal speed of the vehicle is found to always take the maximum value admissible. The actual trajectory is computed via backtracking. This approach allows for both stationary and mobile obstacles, and is able to compute an optimal departure time for the vehicle. The use of generalized gradients and co-states overcomes the hypothesis of regularities of the level set. This promising method is at present still lacking an operational implementation. 

Monte Carlo methods discard exact solutions in favour of faster solutions. Also, they provide a viable technique for fulfilling multiple and competing objectives. A class of Monte Carlo methods makes use of genetic algorithms. They start with guessed routes ("chromosomes") whose subparts ("genes") cross each other and mutate in a random way, in order to find a new route ("offspring") that better fits the objective function of the actual problem. The use of Monte Carlo methods in the context of multi-objective optimization is reviewed in Konak et al. (2006), while an application to ship routing is provided by Szlapczynska (2007). There is also a simulated annealing approach to ship routing (Kosmas and Vlachos, 2012). In this case, in order to find a global optimum a trial route is perturbed in a statistical-mechanical fashion. Given that in Monte Carlo methods there is no exact analytical solution, additional criteria are needed in order to decide whether a solution is satisfactory ("convergence test"). Harries et al. (2003) present an example of a hybrid method making use also of third-party optimization softwares. They employ swell forecasts by ECMWF for the Atlantic Ocean and represent the ship route in terms of parametric curves (B-splines), that are perturbed with respect to the calm sea route. The method relies on the modeFRONTIER package for multi-objective (least time and fuel consumption) optimization. Also, the vessel hydrodynamics are not solved internally, but via the SEAWAY package. Route optimization is claimed just for the open-sea part of the route, and one of their results even shows that the route does not always avoid landmass. 

In discrete methods, the spatial domain is represented by some kind of grid (regular or not) and the optimization is based on recursive schemes. A key concept is the so-called principle of optimality: given a point on the optimal trajectory, the remaining trajectory is optimal for the minimization problem initiated at that point (Luenberger, 1979). This property can be stated as a recursive relation, called "Bellman's condition" in the framework of discrete methods. In Zoppoli (1972) a dynamic programming method for the computation of a least-time ship route in the Indian Ocean is used. The algorithm is able to ingest time-dependent environmental fields by evaluating them at the nearest quantized time value. However, the actual case study provided in the paper just uses stationary fields. Ship operating costs for transatlantic routes are minimized in Chen (1978), where a terminal cost is also included in the objective function. The grid used however is just a band of gridpoints along the rhumbline track, and thus is limited in terms of application when there are complex topological constraints, such as in a coastal environment. Takashima et al. (2009) use dynamic programming for computing minimum fuel routes of a given duration. The propeller revolution number is kept constant during the voyage and its value is adjusted in order to reach the target route duration. The ship course is varied in order to exploit ocean currents. However, the algorithm uses static environmental information, and re-routing is run every 3 h in order to deal with dynamic currents. The dynamic programming method by Wei and Zou (2012) is used to minimize fuel consumption. Both throttle and heading of the vessel can be optimized, again with grid limitations as in Chen (1978). Montes (2005) employs Dijkstra's algorithm to compute least-time routes in time-varying forecast fields. However, the effect of weather on vessel speed is parametrized in terms of subjective parameters ("speed penalty function").
section_name: Review of literature

section_content: There are several recurrent shortcomings in the ship routing literature: the limited capability to deal with complex topological conditions, such as in the coastal environment (Bijlsma, 1975; Hagiwara, 1989; Szlapczynska and Smierzchalski, 2007) ; the need for heuristics or subjective parameters in the optimization algorithm (Kosmas and Vlachos, 2012; Montes, 2005) ; non-explicit use of time-dependent environmental information (Hamilton, 1962; Zoppoli, 1972; Takashima et al., 2009) ; limitations on the functional dependence of the vessel response function (Perakis and Papadakis, 1989; Techy, 2011); and the not yet demonstrated use in an operational environment (Lolla et al., 2014). 

All these issues need to be addressed simultaneously by a model aimed at feeding an operational system that also works in coastal waters, for a wide class of vessels and environmental conditions, taking into account navigation safety according to the latest international standards. In VISIR-I all the above-mentioned shortcomings are overcome. The method is based on an exact graph search algorithm, modified in order to manage time-dependent environmental fields and voluntary vessel speed reduction. It is validated against analytical results. In addition, the graph grid is designed to deal with the topological requirements of coastal navigation. VISIR-I also includes a dedicated motorboat model, and safety constraints for vessel intact stability are considered. 

All these features are described in detail in what follows.
section_name: Our contribution

section_content: In this section we present the method employed by VISIR-I for solving the route optimization problem. First, the problem is formally stated (Sect. 2. 1), then the solution algorithm (Sect. 2. 2), the mechanical model of the ship (Sect. 2. 3) and the processing of the environmental analysis or forecast fields affecting the ship dynamics (Sect. 2. 4) are presented. 

The structure of the computer code is provided in Sect. 2. 5 and a validation of the resulting optimal routes is given in Sect. 2. 6.
section_name: VISIR-I method

section_content: The mathematical problem addressed and solved in an operational way by VISIR-I can be stated as follows. 

A ship route is sought departing from A = (x A, t A ) and arriving at B = (x B, t A + J ) and minimizing the sailing time J defined by
section_name: Statement of the problem

section_content: where x = [x(t), y(t)] T within a set ⊂ R 2 denotes horizontal position, t is the time variable, and

with vessel speed c in calm weather conditions and sustained speed v(x, t) in specific meteo-marine conditions, is the "refractive index" of a horizontal domain of linear extent ds such that ds 2 = dx 2 + dy 2. 

(3)

Note that the integrand in Eq. ( 1 ) can be interpreted as an effective optical depth of the ds wide domain. The notation is reminiscent of the problem of determining the path of light moving in a non-homogenous medium. Indeed light propagates over paths of stationary optical depth (Fermat's principle). Ship speed v results from a dynamic balance between forces and torques acting on and from the vessel. This speed is normally found as the solution of differential equations. However, under steady conditions they reduce to algebraic equations of the type:

where p s is a set of ship parameters and p e is a set of values of relevant environmental fields evaluated at (x, t). Navigational safety also poses limitations on the admissible solutions of Eq. ( 4 ). Such limitations are represented as a set of inequalities of the type:

Parameters p s and p e employed in Eqs. ( 4 ) and ( 5 ) are listed in Table 6. 

If set is also a connected domain, the existence of a solution to the problem stated in Eqs. ( 1 )-( 5 ) entirely depends on Eqs. ( 4 ) and ( 5 ): the quality of the route, specifically its topological and nautical characteristics, is determined by these two equations alone. 

Speed v resulting from Eqs. ( 4 ) and ( 5 ) defines the Lagrangian kinematics of the route:

In order to account for uncertainty in the representation of v, a random noise term could be added to the r. h. s. of Eq. ( 6 ). The problem of finding the least-time route in any meteomarine conditions is thus equivalent to the minimization of J functional with a specified refractive index n(x, t), for assigned boundary values A and B. 

If the time dependence in refractive index n is neglected, the general solution of this problem is known from geometrical optics, with routes being refracted towards optically more transparent regions, according to Snell's law. However, whenever the timescale for changes in the environmental fields is comparable or shorter than the typical route duration, such time dependence can no longer be neglected and new kinematical features of the least-time route may appear. Indeed, it could be advantageous to voluntarily decrease the speed during navigation, as shown in Sects. 2. 2. 2 and 2. 2. 3, or even to wait for some time at the departure location before leaving.
section_name: G. Mannarini et al.: VISIR-I: least-time nautical routes

section_content: The first component of VISIR-I presented here is the shortest-path algorithm. The term "shortest path" is used both in the literature and hereafter with a more general sense than a direct reference to the geometrical distance. Indeed, "shortest" may refer to the spatial or temporal distance, as well as the cost or any other figure of merit of the optimal path.
section_name: Shortest-path algorithm

section_content: Let us consider a directed graph G = [N, E]. In VISIR-I the nodes N are part of a rectangular mesh with constant spacing in natural coordinates (1/60 • of resolution in both latitude and longitude). As shown in Fig. 1, each node is linked to all its first and second neighbours on the grid, forming the set of edges E. Thus, neglecting border effects, there are 24 connections per node. The specific edge arrangement leads to resolve angles of

Whether such 24-connectivity should be increased further is questionable, given that in the present case the environmental analysis and forecast fields are provided on a coarser grid (by about a factor of 4) than the spatial resolution of the graph; see Sect. 2. 4. 

In VISIR-I, the resulting graph is first screened for nodes and edges on the landmass. An edge is considered to be on the landmass if at least one of its nodes is on the landmass or if both nodes are in the sea but the edge linking them intersects the coastline. In such a case, the edge is removed from E, which locally reduces the original 24-connectivity of the graph. When applied to a 1/60 • grid for the Mediterranean Sea region (mode 1 of Fig. 8 ), this procedure still leaves more than 20 million sea edges in E; see Table 2. However, for the actual route computations (mode 2 of Fig. 8 ), just a subset of the whole spatial domain is considered. This subregion is chosen to be large enough so that a further increase in size does not reduce the total sailing time J. At present, the se- The value of the angle θ 12 is provided by Eq. ( 7 ). lection of the subregion shape and extent is left to the user of the model.
section_name: Spatial discretization

section_content: Given that environmental conditions change over a timescale comparable with or shorter than the vessel route duration, edge weights cannot be considered as constants. Thus, in order to solve Eqs. ( 1 )-(3), VISIR-I employs a time-dependent algorithm. 

With reference to the nomenclature in Table 1, a timedependent graph G(t) is fully defined by the sets of nodes, edges, and time-dependent edge weights:

Edge weight a j k ( ) between nodes j and k at time step is defined as

where v j k ( ) is the edge mean ship speed, depending on the average j k of the values of the environmental fields at nodes j and k:

evaluated at time t = t 1 +δ t ( -1). Here t 1 is departure time and δ t is the time resolution of the environmental fields. The functional dependence of v j k ( ) on j k results from the actual model of the vessel, and is derived in Sect. 2. 3. Thus, in VISIR-I, edge weights a j k ( ) are non-negative quantities with a dimension of time ("edge delays") and are time-dependent. Note that Eq. ( 8 ) is the discrete counterpart of Eq. ( 6 ), as long as velocity is non-null. 

There are various methods for computing shortest paths on a graph. For an overview, see Bertsekas (1998) and Bast et al. (2014). A large amount of literature deals with applications for terrestrial networks; see, e. g. Zhan and Noon

Table 2. Parameters of the graph for the Mediterranean Sea after the removal of nodes and edges on the landmass (GSHHG coastline used). In the actual route computations, just a subdomain of the whole basin is selected. Due to border effects, the connectivity ratio A/N < 24 (value that would be expected from Fig. 1 ).
section_name: Time-dependent approach

section_content: Top-left corner latitude 45. 814 Zeng and Church (2009) and Goldberg and Harrelson (2005). 

A key concept in graph methods is the node label, which can be either temporary or permanent. The permanent label X j of node j is the minimum value of the objective function (e. g. J of Eq. 1) attainable at that node. A temporary label Y j is any value before the node label is set to its permanent value. When all node labels are set to their permanent value, Bellman's relation holds (Bertsekas, 1998). 

Depending on the way node labels are updated, graph algorithms may be classified into label setting or label correcting algorithms. A label setting single-source singledestination algorithm with fixed departure time is used here. 

The fact that in VISIR-I destination node is assigned (through x B in Eq. 1) leads to a possible degeneracy of the problem, with multiple shortest paths between the specified source and destination node. In Yen (1971) an algorithm is presented for finding several simple shortest paths. In VISIR-I it is deemed that, in presence of time-dependent environmental fields, it is unlikely that an alternative route with exactly the same navigation time exists. Thus, just the leasttime route is sought after. 

In general, the fact that a graph is time-dependent implies that the shortest path can have special features. In fact, under specific circumstances, the strategy of traversing an edge as soon as possible does not always lead to the shortest path. Also, the shortest path may not be simple (there may be loops) or even not concatenated (Bellman's optimality not fulfilled). This has consequences on the class of algorithm to be applied. Orda and Rom (1990) show that in this respect the critical condition is how fast edge delays vary in time. If a j k (t) is a differentiable function of time t, the authors show that, provided

the best strategy for recovering a shortest path is traversing edge (j k) without waiting at node j (first-in first-out or: FIFO). Indeed, waiting for a time dt > 0 would in best case be compensated but never overcome by a related decrease |da j k | ≤ dt in edge delay. The authors also show that a FIFO time-dependent algorithm has the same computational complexity as a static one. Condition Eq. ( 10 ) may be violated for instance during the decaying phase of a rapidly moving storm. The FIFO condition Eq. ( 10 ) is checked at each run of the model and is generally found to be fulfilled, Fig. 2. Thus, Dijkstra's static algorithm (Dijkstra, 1959) is modified according to the guidelines of Orda and Rom (1990) s FIFO time-dependent algorithm. Related pseudocode is provided in Appendix A. 

Before the algorithm is run, edge delays a j k ( ) are checked for nautical safety constraints, Eq. ( 5 ). If at time step an edge (j k) is unsafe for navigation, we set a j k ( ) = ∞. As seen from Fig. 2, this approach generates gaps in a j k (t) as a function of continuous time t. Such gaps are specific time windows during which the edge is not available for linking its nodes. Whenever edges are removed at specific time steps, a FIFO strategy is no longer guaranteed to be optimal, even though edge delays vary slowly. A source-waiting strategy may be necessary in this case (Orda and Rom, 1990). As a consequence, a route retrieved through a FIFO algorithm may still be sub-optimal. This advanced issue is left open for future versions of the system.
section_name: Parameter Value Units

section_content: As seen above, VISIR-I's strategy regarding navigational safety is to remove unsafe edge delays from the graph by setting their edge weight to ∞, prior to the computation of the optimal route. In addition, as will be shown in Sect. 2. 3. 3, vessel speed v affects the safety constraints. Thus, a modification of v may help in keeping an otherwise unsafe edge in the graph. This, in turn, may contribute to optimization, since avoiding the removal of elements from set A(t) can only lower the length of the shortest path. Such voluntary variations in speed should be contrasted with an involuntary speed reduction due to vessel energy loss, caused by interaction with the environmental fields; see Sect. 2. 3. 2. 

VISIR-I defines, for a vessel with maximum engine power P max, a set of possible values P (s) /P max of engine throttle:

Then, at each edge, speeds v (s) j k ( ) are computed using the ship model. The function g(s) is chosen in order to linearly space engine throttle values; see Table 3 (due to the non-linearity of the vessel model, this choice does not imply linearly spaced values of sustained speed; see Fig. 5 ). Next, throttle-dependent edge weights a (s) j k ( ) are computed via Eq. ( 8 ). Each of these edge weights is checked to see whether it complies with navigational safety constraints. If an edge is unsafe, its edge weight is set to ∞. Finally, the throttle value s * leading to the minimum edge weight is chosen by the algorithm:

and the edge weight is set to such a minimum value:

Given the ordering in Table 3, if s * > 1 then voluntary speed reduction is useful for recovering a faster route which is still safe with respect to ship stability constraints.
section_name: Voluntary speed reduction

section_content: The second component of VISIR-I is a ship model describing vessel interaction with the environment (specified by the forecast fields of Sect. 2. 4) and its stability requirements. 

The following presentation comprises of a balance equation for the propulsion system in Sect. 2. 3. 1, a parametrization of the hull resistance due to calm and rough sea in Sect. 2. 3. 2, and a set of dynamic conditions for the intact stability of the vessel Sect. 2. 3. 3.
section_name: Ship model

section_content: Motorboats are the focus of VISIR-I route optimization. 

For these vessels, propulsion is provided by a thermal engine burning fuel and delivering a torque to the shaft line and, when present, to a gearbox (Fig. 3 ). This torque is eventually transmitted to a propeller, converting it into thrust available to counteract resistance to advancement (Journée, 1976; Triantafyllou and Hover, 2003). 

A full modelling of this energy conversion mechanism is a highly complex task involving, just to mention a few, the efficiency of each of these conversion steps, the effect of hull-generated wake on propeller efficiency and corresponding thrust deduction, and the load conditions of the engine (MANDieselTurbo, 2011). A quantitative description of these processes requires a detailed knowledge of engine, propeller, and hull parameters. This could be obtained by standard measurement procedures such as those provided by the International Towing Tank Conference (ITTC, 2002 (ITTC,, 2011b)). 

For the purposes of VISIR-I, it was deemed sufficient to derive the vessel response function from a balance of thrust and resistance at the propeller. That is, given the brake power P, the total propulsive efficiency η and the total resistance R T applied to the vessel, it is required that

where v is the ship velocity in steady conditions, p s is a set of ship parameters, and p e is a set of relevant environmental field values as in Table 6. One of the possible representations of R T is derived in Sect. 2. 3. 2. Since we are not presently addressing the issue of fuel consumption, the engine rotation speed (rpm) -for which a torque equation is necessary -is not considered. The l. h. s. of Eq. ( 14 ) represents the effective power available at the propeller. The efficiency η results from the product of several components related, for example, to hull shape, propeller, and shaft characteristics (MANDiesel-Turbo, 2011). At the present stage of modelling, the value of η is estimated to a constant (see Table 4 ) and will be refined when a more detailed vessel model is used.
section_name: Propulsion

section_content: In this paper we restrict our attention to displacement vessels. Indeed high-speed planing hulls are characterized by a different dynamic behaviour and deserve a more sophisticated treatment (Savitsky and Brown, 1976). 

When underway, a displacement vessel is subject to various forces hindering its motion. A possible decomposition of the resulting force is to distinguish calm water resistance R c from resistance R aw due to only sea waves,

Each of the addends is meant as the force component opposite the motion of the vessel. The module of the calm water resistance is usually given in terms of a dimensionless drag coefficient C T defined by the equation

where also sea water density ρ and ship's wetted surface S appear. 

As outlined in ITTC (2011a), C T depends not just on viscous effects but also on energy dissipated in gravity waves generated by the vessel ("residual resistance"). The latter introduces a dependence on Froude number Fr which, under Froude's hypothesis, is additive: Fr), where R is Reynold's number and C R is the residual resistance drag coefficient (Newman, 1977). 

For specifying the drag coefficient C T, the statistical method by Holtrop (1984) involves 12 geometrical parameters of the hull. This approach may still imply significant inaccuracies. Indeed, as optimization studies demonstrate, substantial improvements in vessel performances can be achieved through some minor changes to the hull shape, while keeping constant the principal hull parameters (Peri et al., 2001). Hence, it is believed that the most reliable way to account for all the aspects of calm water resistance (both Table 4. Parameters of the ship model. The numerical factor in the formula for Fr value accounts for the conversion of speed v from kt to m s -1 ; g 0 = 9. 80665 is the standard gravitational acceleration; The values of η and ϕ 0 are just guesses. The value of ρ is taken from Cessi et al. (2014) The nautical resistances have the dimension of a force and their unit is the kilo-Newton (kN). frictional and residual) and added resistance in waves would be to use towing tank data for the specific hull geometry, properly transformed to account for scaling effects. However, it is our aim that VISIR-I runs without specifying too many vessels parameters. Thus, C T is taken as a constant. In particular, the C T S product is obtained by equating the maximum available power at the propeller to the power dissipation occurring at top speed c in calm water:
section_name: Resistance

section_content: The impact of assuming a constant C T is to overestimate it at low speeds, as this coefficient is identified using the top speed regime, Eq. ( 17 ). This is quantified in Appendix B, where a sensitivity test is provided, based on a comparison between a constant and a polynomial C T. The contribution of hull fouling to calm water resistance is a long-term timedependent effect and is also neglected. 

In addition to calm water resistance, sea waves are an additional source of ship energy loss (Lloyd, 1998). Various authors have found that wave-added resistance R aw depends on reduced wave number L/λ, where L is ship length. Both radiation (energy dissipated due to heave and pitch movements) and diffraction (energy dissipated by the hull to deflect short incoming waves) contribute to this additional resistance. Both effects were modelled by Gerritsma and Beukelman (1972) in head seas, which however are the most severe conditions in terms of added resistance. They found that diffraction delivers and additional contribution to radiation-induced resistance just for L/λ > 1. In the framework of a comprehensive study of experimental results and several different theoretical methods, Ström-Tejsen et al. (1973) endorsed the method by Gerritsma and Beukelman (1972). However, there is no simple formula which gives the Table 5. Database of vessel propulsion parameters and principal particulars used in this work. See Fig. 3 for the meaning of the geometrical parameters. V1 is a ferryboat while V2 is a fishing vessel. Most data stem from www. marinetraffic. com; T R is estimated from the metacentric height GM using Weiss' method for small roll angles as reported in Benedict et al. (2004) and adding an extra 20 % to account for roll stabilization. Metacentric height is assumed to be GM = 2T /3. is not used by VISIR-I and is provided just for the sake of reference. added resistance in waves for all ship types with good accuracy (Bertram and Couser, 2014).
section_name: Symbol

section_content: In VISIR-I, following the cited literature, a reduced nondimensional resistance σ aw is introduced:

where α is the angle between wave direction and vessel direction of advance (as seen in Fig. 3, α = 0 in case of head waves). The relation between wave amplitude ζ and significant wave height H s is 2ζ = H s. For vessel beam B and draught T see also Table 5. In Eq. ( 18 ) a factor ϕ is highlighted, containing the spectral and angular dependencies. This factor is eventually set to a constant value ϕ 0. This approximation is also done in view of the fact that the full wave spectrum is not used for weighting R aw, as instead done, for example, in Ström-Tejsen et al. (1973). In line with dropping the α dependence in ϕ, the angular dependence of R aw, is ignored by assuming that this force is always opposite to the ship's forward speed in a longitudinal direction (α = 0). Empirical methods are often used for deriving σ aw when the hull geometry is not available in its entirety. They make use of experimental data from a variety of vessels that are fitted in terms of a few parameters, usually the principal particulars. An analysis of the statistical performance of different empirical methods with respect to a database of almost 50 vessels is carried out in Grin (2015). It is distinguished among different L/λ regimes and, where possible, among various ship headings, finding relative errors with respect to the experimental tests in the range of 20-60 %. However, the formulas of these methods are not fully disclosed. Alexandersson (2009), basing on the Gerritsma and Beukelman (1972) method (radiation part only), computes the peak values of the wave added resistance for a database of large ships. He then makes a regression analysis, employing principal particulars 5, respectively. and ship speed. We use its results in a slightly modified way:

Further details of this derivation can be found in Appendix C. The increase of peak value of σ aw with Fr is observed also in the results by Grin (2015), while this is not the case for the increase with L and the decrease with B. Combining Eqs. ( 19 )-( 20 ) with Eq. ( 18 ) shows that an increase in either ship beam or draught leads to an increase in resistance, while an increase in length has the opposite effect. This conclusion should be validated through towing tank measurements on the specific hull geometry. Substituting Eqs. ( 15 )-(20) into Eq. ( 14 ), the following expression is found to relate ship speed to brake power, geometrical vessel parameters, and environmental fields:

where the coefficients are given by

Note that Eq. ( 21 ) is in the form of Eq. ( 4 ) with parameters p s and p e as in Table 6. Sustained speed v is the sole positive root of cubic equation Eq. ( 21 ) (in fact, both k 3 and k 2 coefficients are positive quantities). This root is computed through an analytical expression whose numerical implementation is provided by Flannery et al. (1992, Sect. 5. 6). In Fig. 4 corresponding sustained Froude numbers Fr are displayed. Fr follows a half-bell-shaped curve, with a nearly hyperbolic (∼ 1/H s ) dependence for large significant wave height. While in the results shown by Bowditch (2002, Fig. 3703 ) for a commercial 18-knot vessel, a change of convexity of the Fr curve is not visible (at least for the H s range shown), it is clearly apparent in the results shown by Journée (1976, Figs. 6, 10). 

Our results also prove that, by varying engine throttle, sustained speed does not vary by the same factor at all H s, Fig. 4. This result could not be obtained by factorizing throttle dependence, as in the ship model by Perakis and Papadakis (1989). 

Furthermore, by comparing performances of vessel V1 (ferryboat) and V2 (fishing vessel), it can be seen that the former sustains a larger fraction of its top Froude number at any given significant wave height. This different dynamic behaviour is mainly related to the maximum engine brake power P max of the two vessels. This is found by swapping just P max of the two vessels and keeping the other parameters provided in Table 5 unchanged (not shown). Figure 5 shows how the throttle needs to be adjusted to sustain a given speed in different sea states. An increase in speed requires an overproportional increase in throttle. Lloyd (1998) makes the assumption that the engine delivers constant power at a given throttle setting, regardless of the increased propeller load due to rough weather (note that propeller load is not considered in VISIR-I either). He then finds that the power required for sustaining a given speed steeply rises with wave height (Lloyd, 1998, Fig. 13. 5), in a way similar to Fig. 5. The constantpower hypothesis of Lloyd (1998) is compatible with a turbine engine, which is one of the cases considered in Journée (1976). From Eq. ( 11 ), it follows that in VISIR-I, at constant engine setting (throttle), delivered power is a constant. Thus, it is to be expected that VISIR-I and Lloyd (1998)'s results are qualitatively comparable, as it is indeed found. 

The comparison between V1 and V2 also shows that the two vessels behave quite differently in extreme seas, whereby 4 ) and ( 14 ), and in the inequalities for the safety constraints Eq. ( 5 ). Derived parameters such as T E, σ aw and Fr are omitted. For an explanation of symbols; see (b) (a)

Figure 5. Engine throttle needed for sustaining a given Fr, in calm water resistance only ("H s = 0") and both calm and wave added resistance ("max H s ", i. e. at maximum significant wave height seen in Fig. 4 ). Throttle values correspond to those of Table 3. Panel (a) and (b) refers to ship parameters for vessel V1 and V2 in Table 5, respectively. vessel V1 (the ferryboat) is able to reach more than 30 % while V2 (the fishing vessel) reaches less than 20 % of its top Fr. 

Resistances are evaluated from the sustained speed v as

and corresponding values are shown in Fig. 6. While calm water resistance R c does not explicitly depend on significant wave height H s, R c depends on ship speed which, through Eqs. ( 21 )-( 23 ), depends on H s. Thus, assuming maximum throttle, a functional dependence R c = R c (H s ) can be computed and is displayed in Fig. 6. Due to the fact that k 3 is independent of H s (Eq. 22), calm water resistance R c is dominated by the v = v(H s ) relationship seen in Fig. 4. 

Wave added resistance R aw as a function of H s initially grows quadratically and, for higher waves, only linearly, Fig. 6. This is due to the combined effect of the quadratic dependence on wave amplitude in k 2 (Eq. 23) and the nearly hyperbolic ship speed reduction for large H s seen in Fig. 4. The same trend is observed in (Lloyd, 1998, Fig. 3. 13 ) and Nabergoj and Prpić-Oršić (2007). 

In comparison to V2, vessel V1 exhibits larger resistances. However, for both vessel classes, the R c and R aw curves form "scissors", which are wider for the larger vessel (V1), Fig. 6. This qualitative behaviour compares well to (Journée, 1976; Fig. 12).
section_name: Symbol

section_content: The ship model described so far needs to be complemented by navigational constraints in order to reduce dangerous or unpleasant movements for the ship itself, the crew and cargo. 

Such situations cannot simply be ruled out by designing a vessel in accordance with the Intact Stability (IS) Code, IMO (2008). In fact, specific combinations of meteorological and sea-state parameters may lead to dangerous situations even for ships complying with such mandatory regulations (Umeda, 1999; IMO, 2007). Furthermore, in Belenky et al. (2011) the point is made that new ship forms can make the prescription of the IS code obsolete. This led to the development of "second generation" stability criteria, which is more physics and less statistics based than IS criteria. Computations of this type have recently been carried out by Krueger et al. (2015) for Ro-Ro passenger ships. 

VISIR-I checks for three modes of stability failure: parametric roll, pure loss of stability, and surfriding/broachingto. The theoretical hints below are mainly based on Belenky et al. (2011), while the implementation of the stability checks follows the operational guidance by IMO more closely (IMO, 2007). Because of the limited angular resolu- tion of the graph (Sect. 2. 2. 1), in VISIR-I stability in turning (Biran and Pulido, 2013 ) cannot be taken into consideration, and an unlimited vessel manoeuvrability (IMO, 2002) has to be assumed. 

A realistic assessment of stability failure would require a detailed knowledge of ship hull geometry. In the current version of VISIR-I, however, just principal particulars of the vessel (length, beam, draught) are employed. In addition, even vessel-internal motions and mass displacements, such as the positioning of catch within a fishing vessel (Gudmundsson, 2009 ) and fuel sloshing (Richardson et al., 2005) may have an amplifying effect on the loss of stability. Thus, the bare application of safety constraints described in the following cannot guarantee navigation safety, and the shipmaster should critically evaluate the resulting route computed by VISIR-I, also taking into account the meteo-marine conditions actually met during the voyage and the specific vessel response. While the actual functional form of the safety constraints may be different from what has been implemented, the VISIR-I code addresses the problem of implementing multiple constraints in a numerically efficient way. If necessary, the user can individually switch off such stability constraints by changing the corresponding flags in the namelist file. 

In the following sections, we use the deep water approximation of the wave dispersion relation in order to gain a rapid estimation of the threshold conditions. We can thus estimate the wavelength λ as (T w is the peak wave period) and the wave phase speed or celerity c p as

Then, assuming a fully developed sea (Pierson-Moskowitz spectrum), the wave steepness can be estimated as

This result can be inferred from the plot of characteristic seas reported by Ström-Tejsen et al. (1973). Wave steepness is larger than the value obtained in Eq. ( 28 ) for partially developed seas and smaller for dying seas.
section_name: Stability

section_content: When a ship is sailing in waves, the extent of the submerged part of the hull changes in time. For most hull shapes, this also involves a change in the waterplane area. This in turn influences the curve for the righting lever (GZ), which is fundamental to ship stability. Indeed, if wavelength λ is comparable to ship length L and waves are met at a specific frequency, the change in GZ may trigger a resonance mechanism, leading to a dramatic amplification of roll motion (Belenky et al., 2011). A famous naval casualty ascribed to this mechanism of stability loss is reported in France et al. (2003). 

The mathematical formulation of parametric roll is based on the solution of Mathieu's equations and the computation of Ince-Strutt's diagram. It shows that parametric roll occurs when encounter wave period T E satisfies the condition

where T R is the ship's natural roll period (Spyrou, 2005) and the ± sign in Eq. ( 29 ) accounts for both head and following seas. 

In VISIR-I the encounter period T E is obtained by applying a Doppler's shift to T w and reads

where Fenton's factor K defined by Eq. ( 47 ) is used and v is given in knots. Instead, IMO's formula for T E provided in IMO ( 2007 ) corresponds to the deep water approximation, i. e. to the case K = 1. Since in shallow waters and large wave periods K < 1, IMO's formula may lead to an overestimation of T E. Levadou and Gaillarde (2003) observe that a smaller GM also implies a larger natural roll period T R and thus a parametric roll experienced in presence of longer waves. Spyrou (2005) points out that, while any encounter angle α can in principle lead to parametric roll, vessels with low metacentric height GM (and thus large T R ) may be more prone to experience parametric roll during following than head seas (due to larger |T E |). 

Following Levadou and Gaillarde (2003) and the wave height criterion reported for L < 100 m in Belenky et al. (2011), the parametric roll hazard condition is implemented in VISIR-I as

together with Eq. ( 29 ) expressed in the form of the following inequalities:

where the coefficients in Eqs. ( 33 )-( 34 ) should be related to the roll damping characteristics of the vessel (Francescutto and Contento, 1999), but for the current version of VISIR-I they are taken from Benedict et al. (2006). Formula Eq. (30) shows that T E period may be tuned by varying the speed and course of the vessel. Thus, to prevent parametric rolling, a routing algorithm may suggest either a voluntary speed reduction or a route diversion. As shown in Sect. 2. 2. 3 and as will be seen in the case studies (Sect. 3), VISIR-I is able to exploit either option.
section_name: Parametric roll

section_content: This mode of stability failure is triggered by a similar condition to the parametric roll. However, it does not involve any resonance mechanism and thus may be activated by a single wave. In fact, if the crest of a large wave is near the midship section, stability may be significantly decreased. If this condition lasts long enough (such as during following waves and a ship speed close to wave celerity), the ship may develop a large heel angle, or even capsize. 

According to Belenky et al. (2011) a useful criterion for distinguishing ships prone to pure loss of stability involves a detailed knowledge of hull geometry. The IMO guidance (IMO, 2007), however, suggests using just ship-wave kinematics. This is also the criterion adopted in VISIR-I and can be stated as the following conditions to be simultaneously verified:

where ship speed v is given in kt. 

Using also Eqs. ( 26 )-( 27 ) it can be seen that Eq. ( 38 ) implies (for exactly following seas) a sustained speed v between 43 and 67 % of wave celerity c p.
section_name: Pure loss of stability

section_content: Surfriding is the condition where the wave profile does not vary relative to the ship. That is, the ship moves with a speed equal to wave celerity: v = c p. In this case, the ship is directionally unstable, with the possibility of a sudden and uncontrollable turn known as "broaching-to". 

The simplest modelling of this mode of stability failure starts with the computation of the force of the wave-induced surge which is able to balance the difference between total resistance and thrust provided by the ship. A critical point may then be reached, where surging is no longer possible and the ship is captured by the surfriding mode (Belenky et al., 2011). This phase transition is a heteroclinic bifurcation (Umeda, 1999). 

In IMO (2007) a surfriding condition is proposed which just takes into account ship speed and length, independently of wave steepness. Based on numerical simulations, Belenky et al. (2011) overcomes this simplification, with the finding that the phase transition is less likely for less steep waves. 

In VISIR-I, the following surfriding hazard criteria reported in Belenky et al. (2011) are considered:

where the critical Froude number is given by

Using Eq. ( 28 ) its typical value is found to be Fr crit = 0. 31. Condition Eq. ( 40 ) was added to VISIR-I since Fr crit is reported in Belenky et al. (2011) just for the range Fr ∈ [1/40, 1/8]. Condition Eq. ( 42 ) was complemented with an α dependence in analogy with Eq. ( 38 ) in order to account for following-quartering seas. This implies that surfriding is less likely to occur for quartering than following seas, since Fr is multiplied by a factor which may be as small as 1/ √ 2. 

Of note is that all VISIR-I safety constraints described above, Eqs. ( 31 )-( 42 ), are implemented in negative, i. e. as the set of conditions possibly leading to a stability loss. Nevertheless, they are all still in the form of Eq. ( 5 ) with parameters p s and p e as in Table 6.
section_name: Surfriding/broaching-to

section_content: We distinguish the environmental fields between static (bathymetry and coastline) and dynamic fields (waves, winds, currents). In VISIR-I, bathymetry and coastline are employed to ensure that navigation occurs in not too shallow waters and far from obstructions. Of the dynamic fields, just wave forecast fields are used, as explained in Sect. 2. 4. 2.
section_name: Environmental fields

section_name: Static fields

section_content: A 1/60 • (= 1 nautical mile or 1 NM) bathymetry is employed in VISIR-I. The data set (NOAA Digital Bathymetric Data Base 4 ) is used for a twofold purpose: i. Along with the coastline database, bathymetry is needed for computing a land-sea mask for safe navigation. The first step is to select edges (j k) satisfying the condition that edge averaged sea depth z = (z j + z k )/2 is larger than ship draught T :

In other words, just a strictly positive under keel clearance UKC = z -T is admitted for navigation. 

ii. Bathymetry is needed also for a more accurate estimation of wavelength λ, which is an important quantity for vessel stability checks of Sect. 2. 3. 3. Indeed deep water approximation tends to overestimate λ in shallow waters. Instead, VISIR-I employs Fenton's approximation (Fenton and McKee, 1990) which, upon the introduction of the deep water limit λ 0 for the wavelength of the spectrum component of peak period T w,

can be rewritten as follows:

As seen from Eq. ( 47 ), in order for λ to sense the effect of shallow water, z should be small with respect to the scale set by λ 0.
section_name: Bathymetry

section_content: The coastline database is used in VISIR-I for a preliminary removal of graph edges on the landmass (Sect. 2. 2. 1) and, jointly with the bathymetry, for the computation of a nautically safe land-sea mask (see below). 

To this end, the NOAA Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG 5 ) is employed. Just two hierarchical levels are considered: the coastline of the Mediterranean basin and its islands. The minimum distance between coastline data points is variable and is in some cases below 100 m. 

A joint depth-coast land-sea mask is obtained by multiplying the mask defined by Eq. ( 44 ) with a mask of offshore grid points. This way, VISIR is suited for complex topology: domains with presence of peninsulas, islands, and archipelagic seas can all be successfully addressed (see also case studies in Sect. 3). 

Due to the quite different spatial resolution of the coastline and the environmental fields, a regridding procedure is employed for reconstructing the coastal fields. 

1. Fields are extrapolated inshore by replacing missing values of sea fields with the average of the first neighbouring grid points, Fig. 7. Such "sea-over-land" procedure can be iterated in order to define field values on further neighbouring land grid points. This approach is distinguished by the extrapolation used in De Dominicis et al. ( 2013 ) by the number of neighbours used (8 and not just 4) and the absence of the condition that at least two neighbouring grid points have assigned field values. Yet an another procedure is used by Kara et al. (2007) for correcting atmospheric fields from land contamination: a weighted sum over the compact nine-point stencil is computed and the target point is filled if the weights sum to at least a minimum score. 

2. The fields are bi-linearly interpolated to the target grid. 

In VISIR-I this is the bathymetry grid. Thus, spatial resolution of wave fields is enhanced from the original 1/16 to 1/60 •.
section_name: Coastline

section_content: The dynamic environmental fields are used in VISIR-I for the computation of sustained ship speeds and safety constraints. 

In the present version, just the effect of waves is considered, which is deemed to be the most relevant for medium-and small-size vessels. The effect of wind and sea currents is planned for future development. In fact, 1. wind drag may be significant for vessels with a large freeboard and/or superstructure area (Hackett et al., 2006) ; Target grid for the interpolation performed after application of sea-over-land is drawn as a dashed grid (for ease of presentation, it is drawn exactly 4 times finer than the original grid). Also shown in (b) is the land-sea mask of the target grid (green area). 

2. sea current drift is relevant especially in proximity to strong ocean currents (Takashima et al., 2009) and for vessels with large draughts that are not too fast;

3. wave effects include both drift and involuntary speed reduction. The drift is due to nonlinear mass transport in waves (Stokes' drift, Newman, 1977). It is small when the reduced wave number L/λ is smaller than unity and increases significantly when L/λ ≈ 1 (Hackett et al., 2006). Involuntary speed reduction in waves was instead detailed in Sect. 2. 3. 

Thus, the effect of wind drag may be neglected for nottoo-large vessels, and the effect of current and wave drift may be neglected for vessels able to sustain significantly larger speeds than the current magnitude. In addition, since coastal wave fields may be affected by the extrapolation/interpolation procedure, and due to the current resolution of the bathymetry grid (1 NM) (Sect. 2. 4. 1), very small vessels sailing coastwise on short routes should be removed from the scope of this system. Thus, we roughly estimate the range of admissible vessel lengths L to be between 10 and a few tens of metres. 

The current version of VISIR-I employs wave forecast fields from an operational implementation of the Wave Watch III (WW3) model (Tolman, 2009) in the Mediterranean Sea, delivered by INGV (Istituto Nazionale di Geofisica e Vulcanologia) as a part of the Mediterranean Ocean Forecasting System (MFS) system. WW3 is a spectral model that considers (for deep water conditions) as action source and sink terms: wind forcing, whitecapping dissipation, and nonlinear resonant wave-wave interactions. Details on the physical mechanisms implemented in the current application in the Mediterranean Sea can be found in Clementi et al. (2013). Service: Pinardi and Coppini (2010), Oddo et al. (2014 ), and Tonani et al. (2014, 2015). The coupling involves an hourly exchange of sea surface temperature, sea surface currents, and wind drag coefficients between the two models (Clementi et al., 2013). The WW3 model is horizontally discretized on a 1/16 • mesh. Wind forcing is through 1/4 • resolution 6 ECMWF model forecast fields with 3-hourly resolution for the first 3 days and then a 6-hourly resolution. For the case studies of Sect. 3, fields from WW3 run in hindcast mode are employed: ECMWF analyses are used as a forcing for both the wave and the hydrodynamic model and NEMO is run in data assimilation mode. The spectral discretization of the current WW3 implementation is 24 equally distributed angular bins (i. e. 15 • ) and 30 frequency bins ranging from 0. 05 Hz (corresponding to a period of 20 s) to 0. 79 Hz (corresponding to a period of about 1. 25 s). The operational product used as input by VISIR-I, however, does not contain the full spectral dependence, but just the peak wave period T w, significant wave height H s and wave direction θ w. 

Hourly output fields of the MFS-WW3 model are employed by VISIR-I.
section_name: Dynamic fields

section_content: Here we present the main steps in the computational implementation of VISIR-I into a computer code. The code itself and a data sample can be obtained following the instructions provided in Sect. 5. The flow chart in Fig. 8 shows that there are two distinct VISIR-I functioning modes. In both modes, the first step is to prepare the model grid hosting graph nodes and edges. 

Mode 1 is needed to produce the database of nodes and edges neither lying on the landmass nor crossing it; see Table 2. Sea nodes are computed first, since sea edges are a subset of the edges linking sea nodes (an edge can link sea nodes and still cross the landmass). This selection is a timeconsuming process and at the same time completely independent of the forecast fields. Thus, mode 1 is run once for a given topology of the domain (coastline) and graph structure (grid resolution and connectivity). The resulting database of nodes and edges is then employed as VISIR-I runs in mode 2. 

Mode 2 is the functioning mode for the operational use of VISIR-I. First of all, the ship model is evaluated. Equation ( 21 ) is solved and a look-up table of ship speed values v = v(P (s), H s ) as a function of engine power settings P (s) and significant wave heights H s is prepared, as described in Sect. 2. 3. 2. All environmental fields are then restricted to the domain where the route is to be searched. Gridded fields are converted to edge average quantities through Eq. ( 9 ). In order to compute the time-dependent edge weights a j k ( ), the look-up table v = v(P (s), H s ) is linearly interpolated for the actual H s value relative to each edge. At the same time, edge weights of set A(t) that at specific times t are not compliant with the navigational safety constraints are set to ∞. The shortest-path algorithm is then run twice. First, it is run in its time-independent version using the geodetic distance between nodes as edge weight 7 :

This computes a still-safe geodetic route from a topological viewpoint (coastline and bathymetry already checked at previous steps). The time-dependent shortest-path algorithm is then run with time-dependent edge weights a j k ( ) from Eq. ( 8 ). The output of the shortest-path algorithm is a set of nodes and times at which they are visited. This information is necessary and sufficient for reconstructing all environmental fields (H s, θ w, T w, T E, z) and ship status variables (x, P, v, v) along the route. 

In VISIR-I, for long routes, the computing time is dominated by the preparation of the edge weights and the shortestpath computation. The computing time τ for the various model components can be represented by polynomial fits in y ) is shown as a colour background field (the feature at about 10 • 45 E is due to the projection into spherical coordinates of a speed field constant along the x linear coordinate). The geodetic and the optimal routes are displayed with squared black and circular red markers respectively. The analytical solution Eq. ( 50 ) (inverted cycloid) is plotted as a dashed black line. The length and duration of these three routes are compared in Table 10. 

terms of the number N of gridpoints included in the selected spatial domain for the route computation:

with coefficients, as in Table 9. The performance of the shortest-path component of the model could be improved by making use of data structures such as binary heaps (Bertsekas, 1998).
section_name: Outline of the computational implementation

section_content: An exact validation of the optimization algorithm of VISIR-I and the forthcoming post-processing phase is possible in the case of time-invariant fields. However, algorithmic complexity and pseudocode do not substantially differ for the case of time-invariant and time-dependent fields, as pointed out in Sect. 2. 2. 2. In fact, they basically differ just in using edge weight a j k ( ) instead of a j k in row #9 of pseudocode in Appendix A. Thus, a validation of the algorithm for timeinvariant fields covers a more general scope. We exploit the cycloidal curve, being the solution to problem Eqs. ( 1 )-(3) if speed v is proportional to the square root of one of the horizontal coordinates. If speed is given by v = √ 2g(2Ry) the solution is an inverted cycloid:

where 2R is the distance between departure and arrival point along y direction and 0 ≤ y ≤ 2R (Lawrence, 1972). Thus, Table 9. VISIR model performance metrics. The coefficients are identified by least-square fits of Eq. ( 49 ). Routing jobs were run for grids of size N in the range [8×10 1 -8×10 4 ]. The c 0 offsets are constrained for all but the case of the "job total" computing time. Computer features: 3. 5 GHz Intel Core i7 processor with 32 GB RAM, 1600 MHz DDR3. 

Route type c 0

Geodetic 0. 0 1. 0 × 10 -4 3. 2 × 10 -9 99. 9 0. 1 Optimal 0. 0 1. 1 × 10 -4 3. 3 × 10 -9 99. 9 0. 1 Edges 0. 0 9. 6 × 10 -4 7. 7 × 10 -9 98. 9 2. 1 Job total 4. 0 4. 2 × 10 -4 2. 2 × 10 -8 99. 4 4. 9

the aspect ratio of the cycloid is defined solely by parameter R. On the other hand, time J for moving between the two endpoints of the curve under the influence of a "gravity force" also depends on g parameter (see formulas in Table 10). Figure 9 proves that the VISIR-I optimal route follows the analytical trajectory of Eq. ( 50 ). The geometrical differences are due to the connectivity of the graph, leading to the angular resolution θ 12 given by Eq. ( 7 ). In Mannarini et al. (2013), the effect of graph connectivity on the representation of analytical routes was quantified in the absence of environmental fields. Also in terms of navigation time, the VISIR-I optimal route is quite accurate with respect to the cycloid (see Table 10 ). While the length error is about 2 %, the error in route duration is just 1 %. This is because larger misfits with respect to the cycloidal route are found in the lower latitude portion of the route, where the advance speed is highest and thus a relatively shorter time is spent, and this leads to a smaller accumulation of temporal errors. 

Note that the cycloidal profile is compatible with Snell's law of refraction, as the route is refracted in order to reach the optically more transparent (higher speed) region the soonest. Instead, the rhumb line connecting departure and arrival points does not sufficiently exploit such a high-speed region and lags behind by more than 18 %; see Table 10.
section_name: Validation

section_content: In VISIR-I, the choice of the vessel parameters (Table 5 ), the variety of possible sea states, and the freedom to select departure and arrival from any two points in the Mediterranean Sea give rise to a considerable number of route features. In this section we generate a few prototypical situations, demonstrating the features of the model presented in Sect. 2. As mentioned in Sect. 2. 4. 2, analysis rather than forecast fields are used for computing the results shown in this section. The FIFO condition of Eq. ( 10 ) is checked at each time step of the analysis fields. It turns out that FIFO is always fulfilled in the cases considered. 11, in this case, the geodetic route takes longer than the optimal route to reach the destination (d). Animation of the route is provided at http://dx. doi. org/10. 5446/18087. 

Table 10. Cycloidal benchmark: length and duration of the three routes shown in Fig. 9. The free parameters in Eq. ( 50 ) are R = 14. 6 NM and g = 10 -3 m s -2. Quantities in rows 2 and 5 refer to the numerical results using the VISIR-I grid. Quantities in rows 3 and 6 are errors computed with respect to the values in the "Perfect cycloid" column of rows 2 and 5. 5. 

island of Favignana (37. 9 • N, 12. 3 • E) and diverts further southwards while crossing the Strait of Sicily. Finally, after a course change towards starboard, it reaches Tunis. This occurs at a time when an identical vessel on the geodetic route with the same departure time has not yet reached its destination (Fig. 10d ). 

Considering the motion of the wave height field as well, the optimal route attempts to maximize the time spent in calmer seas, where, due to the smaller added wave resistance (Eq. 18), the sustained speed is higher. This is why, though longer in terms of sailed miles, the optimal route is significantly faster than the geodetic route (Table 11 ). 

Figure 11 further analyses the temporal evolution of the two routes. Beginning about t = 6 h after departure, a sawtooth feature in the time history of H s and v variables is displayed. This is due to the temporal variation of the wave field, which is fast on the scale of the time step duration (δ t = 1 h) at which the analysis fields are provided. However, we checked that the FIFO condition of Eq. ( 10 ) is still satisfied for this route. Superimposed on the saw-tooth, there are smaller steps in both H s and v time series, roughly any δ g = δx/v ∼ 0. 1 h, where δx = 1 NM is the VISIR-I graph grid spacing and v ≈ 10 kt is the ship speed at t = 6 h. These smaller steps are due to the strong spatial gradients of the local significant wave height field. The encounter wave Table 11. Summary metrics for the case study routes displayed in Figs. 10-15. Values in bold between brackets for case study #3 refer to the optimal route without voluntary speed reduction (engine throttle forced to always be 100 %). Variations are computed as = 100•(opt/gdt-1). N is the number of sea nodes found in the bounding box selected for the computations; T O is the time spent in the computation of the optimal route; while T J is the total job computing time (excluding rendering of maps and time series). T J is shorter in case voluntary speed reduction is not applied, since edge weights have to be evaluated at just a single engine throttle. T O and T J refer to the performance achieved on a 3. 5 GHz Intel Core i7 processor with 32 GB RAM memory, 1600 MHz DDR3. 

Case 5. 

period panel shows that, at about t = 12 h, T E of the optimal route nearly matches T R. This is one of the necessary conditions for parametric rolling, as required by Eq. ( 34 ). However, the panel with the danger indices shows that such danger condition is not activated. This is due to a large wavelength λ > 2L, not matching criterion Eq. (31).
section_name: Mediterranean Sea case studies

section_content: In the second case study, a transfer of fishing vessel V2 of Table 5 between the islands of Crete and Rhodes (Greece) is assumed to occur during a Meltemi (north wind) situation, typical for the Aegean Sea. 

In Fig. 12 the geodetic and optimal routes are displayed on top of the wavelength field. In this case, wavelength λ is often comparable to vessel length, as clearly seen from the λ time history in Fig. 13. This condition favours, along the geodetic route, the infringement of the stability criteria for both parametric roll and pure loss of stability of Sect. 2. 3. 3. In fact, the reduced wavelength λ/L controls the activation of all safety constraints (see Eqs. 31, 35, and 39). However, this is not the reason for the diversion south of island of Karpathos (35. 4 • N, 27. 2 • E) suggested by VISIR-I, which is still driven by the refraction effect through calmer seas (this can be seen by switching off the safety constraint checks). Both routes correctly avoid all obstructions, maintaining a positive UKC,

(e) (f)

Case study #3. Geodetic (black markers) and optimal (red markers) route from Gibraltar (UK) to Ben Abdelmalek Ramdan (Algeria) for vessel V2 of as required by Eq. ( 44 ). Finally, we note that for such an "archipelagic" domain and large T w, shallow waters significantly affect wave dispersion. For example, the wavelength λ during the first hour of navigation would be overestimated by about 20 % if, in place of Fenton's approximation Eqs. ( 46 )-( 47 ), the deep water approximation were employed in a region where z < 50 m (not shown).
section_name: Case study #2

section_content: In the third case study, a voyage of fishing vessel V2 of Table 5 from Gibraltar (UK) to Sidi Ali (Algeria) is simulated during a wave event past the strait of Gibraltar into the Mediterranean Sea. 

Figure 14 shows a northbound diversion of the optimal route compared to the geodetic route, being instead along a line of constant latitude. The diversion results in the optimal route reaching the destination significantly later than the geodetic route; see Table 11. However, this is still compliant with the least-time objective of VISIR-I, as in this case the northbound diversion is forced by the safety checks (as proven from Fig. 14e, where they are disabled). Both pure loss of stability and surfriding may occur along the geodetic route, as seen from Fig. 15. This is due to prevailing following seas and relatively short wavelength, compared to vessel length. Parametric roll instead is inhibited due to large |T E |, resulting from fetch, Fig. 15b. Also, the relatively small significant wave height leads to a sufficiently high sustained speed for the threshold condition on the Froude number Eq. ( 42 ) to be overcome, thereby originating surfriding conditions for part of the geodetic route, Fig. 15f. This is one of the reasons why for this route the voluntary speed reduction of Eq. ( 12 ) is also at work. As seen from Fig. 15e, the algorithm suggests reducing the throttle to 85 or 70 % (i. e. s = 2 or s = 3 of Table 3 ) for a total of several hours, starting from t ≈ 13 h. This reduces the sustained speed (Fig. 15c ), enabling the vessel to sail with following seas (145 ≤ |α| ≤ 174) without being exposed to surfriding. According to Fig. 15e, engine throttle is reduced and restored again six times in the course of about 10 h. The resulting effect on fuel consumption and onboard comfort is neglected by VISIR-I, as the sole optimization objective is the total time of navigation. Indeed, as seen from Table 11, throttle reduction results in a 20 min faster route than with throttle always kept at 100 %. In the latter case, a southbound diversion is also needed in the last part of the route, as seen from Fig. 14f. As can be seen in Fig. 11, fields along the route ex- 5. 

hibit a saw-tooth feature after t = 10 h, due to the rapidly evolving wave field. However, the FIFO condition of Eq. ( 10 ) is still satisfied. Finally, the impact on this case study of the approximation of constant drag coefficient C T (Sect. 2. 3. 2) is explored in Appendix B.
section_name: Case study #3

section_content: In this paper, we have presented the scientific basis of VISIR-I, a ship routing system, as well as results of its computation of optimal routes in the Mediterranean Sea. The system is designed for flexible modelling of the vessel and its interaction with the environment. Time-dependent analysis and forecast fields from oceanographic models are employed in input. 

The optimal routes computed by VISIR-I were shown to correctly avoid islands and waters shallower than ship draught. Vessel course is generally refracted towards regions of larger sustained speed, allowing in some cases to sail a longer path and reach the destination earlier than along the rhumb line. VISIR-I optimal routes are checked for vessel intact stability, in terms of compliance with IMO regulations and more advanced research results. In some cases, it is these safety criteria, and not the refraction, being responsible for route diversions. The algorithm is also able to compute voluntary speed reductions. The vessel parameters needed to run the model are limited to basic propulsion data and hull principal particulars, making the system accessible for on-demand computations even by non-professionals of navigation. Several issues require further improvements. In relation to the time-dependent algorithm, in the case of rapid changes in the analysis or forecast fields, the optimality of the route retrieved by the model is no longer guaranteed (see FIFO condition in Sect. 2. 2. 2). This however is not the case in all the case studies presented in this paper. In general, the optimal departure time may include a waiting time, which is not handled by the current scheme. Furthermore, the discretization of the dynamic fields may lead to oscillations in the time history of optimal speed, as seen in Figs. 11 and 15. In some cases, the limited angular resolution of the model and the form of the safety constraints may also lead to sudden course changes in the optimal route. For more affordable computations, the graph grid may need a redesigning, thus reducing the density of gridpoints in open seas through the use of a nonuniform mesh. An unstructured (Shewchuk, 2002) or adaptive refinement mesh (Berger and Colella, 1989 ) could be considered. 

For the environmental fields, the most urgent upgrade seems to be accounting for wind, especially for larger ves-sels, and for currents, at least for slower boats. Wind was recently added to a variant of VISIR-I for sailboat routing (Mannarini et al., 2015). Coastal models and limited area weather models could provide the high spatial resolution required for more realistic computations. 

Concerning the ship model, a more advanced parametrization of both calm water and wave added resistance (Sect. 2. 3) could be devised, e. g. by employing data measured in a towing tank. The dependence of vessel drag coefficient on the Froude number, as well as the angular nature and spectral dependence of the total resistance could be considered. Engine modelling could be refined by introducing engine torque and the number of revolutions as well as propeller parameters. A more realistic ship model would also allow for more sensitive optimization objectives, such as onboard comfort. 

A summary of the main approximations employed in VISIR-I is provided in Table 7. 

The first operational implementation of the system took place in the Mediterranean Sea (Mannarini et al., 2016) and can be accessed from the web page http://www. visir-nav. com and the mobile apps linked there. Extension of VISIR-I to any other marine domain is possible. To this end, the corresponding databases for shoreline and bathymetry, along with the forecast or analysis fields are required. Depending on the extension and topological features of the domain, the graph grid and its connectivity deserves a redesign. Furthermore, other environmental fields (such as sea currents, winds, tropical cyclones, sea ice) may also be relevant, depending on geographical domain and vessel class addressed, requiring a revision of the analysis done in Sect. 2. 4. 2 of this manuscript and, correspondingly, an update of the vessel model. 

In the future, VISIR could be generalized to other optimization objectives, such as bunker savings, by suitably modifying the refractive index in Eq. ( 2 ) and adding a torque balance equation in the vessel model in Sect. 2. 3. Another interesting upgrade could be to account for the stochastic nature of the environmental fields. For the vessel modelling, an extension to planing hulls is possible (Savitsky and Brown, 1976). 

In conclusion, we would like to stress the potentiality of VISIR to offer the scientific and technical communities an open platform whereby various ideas and methods for ship route optimization can be shared, tested, and compared to each other. In this respect, the fact that in VISIR-I -through this paper and related source code -the various system components (vessel model, shortest-path algorithm, and processing of the environmental fields) are openly documented and made publicly available should enable unprecedented developments in the efficiency and safety of navigation.
section_name: Conclusions

section_content: The VISIR-I code is made available under the GNU General Public License (Version 3, 29 June 2007) at www. visir-model. net. The VISIR-I code is written in Matlab and can be run on any workstation or laptop. The currently supported architecture is * nix (tested on Mac OS-X 10. 9+ and Linux CentOS 6. 2). Required third party software are the MEXCDF libraries (for reading netcdf analysis/forecast files) and the m_map package (for visualization of maps).
section_name: Code and data availability

section_content: The pseudocode for the time-dependent shortest-path algorithm employed in VISIR-I (see Sect. 2. 2. 2) is provided in this Appendix. 

It is organized into three main parts: initialization of node labels and indices (rows 2-6); main iteration loop (rows 7-10 and 13-17); exit condition (rows 11-12). 

The input arguments are the start and end nodes j s and j e, the set of graph edges {(j k)}, and the set of time-dependent edge weights {a j k ( )}. Index f k in rows 3, 4, 10 represents the predecessor ("father") of node k along the shortest path. δ t is the time step of the time-dependent edge weights. The temporary and the permanent label of node j are respectively Y j and X j. Row 6 implements the FIFO hypothesis by requiring edge weights to be evaluated at the first available time step of the environmental fields, as explained in Sect. 2. 2. 2. In order to speed up access to the set of neighbours of a given node in row 8, a forward star representation of the graph is employed (Ahuja et al., 1988). The minimum label search in row 15 is a typical feature of any Dijkstra's method. The list V contains all nodes k whose permanent label X k is still unset. The Supplement related to this article is available online at doi:10. 5194/gmd-9-1597-2016-supplement.
section_name: Appendix A: Pseudocode for the time-dependent graph method

section_content: Acknowledgements. Funding through TESSA ( PON01_02823 ) and IONIO (subsidy contract no. I1. 22. 05 ) projects is gratefully acknowledged. Mannarini and Pinardi were partially funded by the AtlantOS project ( EC H-2020 grant agreement no. 633211 ). 

Edited by: R. Marsh

section_name: Appendix B: Beyond a constant drag coefficient C T

section_content: In order to numerically evaluate the impact of a constant drag coefficient C T (see Sect. 2. 3. 2), VISIR-I routine ship_resistance. m can be used to solve Eq. 14 in the 25 presence of any polynomial form of C T = C T (v). In particular, we have tested

for various values of q. If the value of γ q is identified at the top powering conditions and H s = 0, (cp. Eq. 17), it reads: 30

where k 3 is given by Eq. 22 and c is the vessel's top speed.
section_name: Appendix B: Beyond a constant drag coefficient C T

section_content: In order to numerically evaluate the impact of a constant drag coefficient C T (see Sect. 2. 3. 2), VISIR-I routine ship_resistance. m can be used to solve Eq. ( 14 ) in the presence of any polynomial form of C T = C T (v). In particular, we have tested

for various values of q. If the value of γ q is identified at the top powering conditions and H s = 0, (cf. Eq. 17), it reads

where k 3 is given by Eq. ( 22 ) and c is the vessel's top speed. 

The ρS dependence is cancelled in the resistance R c :

generalizing Eq. ( 24 ). 

The case q = 0 corresponds to the results shown in Sect. 2. 3. 2, while q = 0 leads to a polynomial of degree q +2 for the residual resistance (we are still neglecting the v dependence of the frictional component in R c ). In the following, we augment the results already provided for q = 0 in Sect. 2. 3. 2, and report a comparison of the q = 1, 2, 3 cases in Fig. 16 and Table B1. 

First of all, we note that, at maximum engine throttle, the vessel speed curve as a function of H s (Fig. 16a, b ) is scarcely affected by the value of q. This is due to the fact that, for H s = 0 and maximum throttle, the speed is constrained to always be c per construction while, for large H s, the wave added resistance dominates the calm water resistance (cf. Fig. 6 ) and consequently the residual resistance. The effect of varying the value of q is also displayed by the plots of engine throttle needed for sustaining a given F r (Fig. 16c, d ). As expected, for calm sea the minimum sustained speed increases with q, since a lower C T -keeping all other parameters fixed -implies a higher vessel speed. 

Finally, we can visualize the effect of q on the route kinematics for the case study #3 of Sect. 3 from panels Fig. 16e, f. Such a case study is chosen for display since the changes to route geometry due to q = 0 are most noticeable. However the other cases were also addressed by the sensitivity test and the results are summarized in Table B1. There is an effect on the length of the diversion of the optimal with respect to the geodetic route. The overall kinematics of the route are also affected, as the same sea state is experienced at (slightly) different times during navigation. From Table B1 it is seen that the total navigation time is reduced for larger q, as expected. Maximum time-savings sum up (for q = 2) to about 7 % of the duration of the q = 0 route (case study #1). Thus, we can conclude that -though a polynomial behaviour of C T will shorten the duration of the routes and could be considered B1. Summary metrics for the routes of all case studies in Sect. 3 and different values of q parameter in Eq. (B1). Voluntary speed reduction is allowed. For both the geodetic and the optimal route, = J (q)/J (0)-1 is the relative difference in navigation time with respect to the case of constant C T (i. e. q = 0). 

Case study Quantity Units Geodetic route Optimal route for the next version of VISIR -the initial approximation of a constant C T does not lead to dramatically different results.
section_name: Appendix B: Beyond a constant drag coefficient C T

section_content: The steps leading to the expressions Eqs. ( 19 )-( 20 ) for the non-dimensional added wave resistance σ aw are described in the following. We start from (Alexandersson, 2009, Eq. 7. 11) results. They are based on Gaussian fits of the outcomes of the method by Gerritsma and Beukelman (1972). A set of seven vessels of different types (Reefer, container, RoRo, tanker) in the range 130-280 m is considered. (In VISIR-I we are addressing smaller vessels, thus such parametrization should be updated in the next versions). The factor of (Alexandersson, 2009, Eq. 7. 11) containing the prismatic coefficient is neglected, since its exponent is very close to zero. The normalized pitch radius of gyration is set to k yy = r yy /L = 1/4. Furthermore, we believe that the longitudinal position of the centre of gravity L CG appearing in Alexandersson (2009, Eq. 7. 11 ) should be replaced with its normalized counterpart L CG /L, 9 that we set to 1/2. The original power law dependence σ aw ∼ Fr 0. 64 is replaced by a linear dependence, Eq. ( 19 ). This is done to retrieve an analytical solution for ship speed v, (see Eq. 21). In order to obtain such linear dependence, a least square fit of σ aw ∼ Fr 0. 64 is forced through the origin over the domain Fr ∈ [0, c/ g 0 L] (C1)

of the independent variable. The range Eq. (C1) of Fr represents the operational regime of the vessel. The slope 1/ Fr of the fitted function identifies the reference Froude number Fr used in Eq. ( 19 ).
section_name: Appendix C: Non-dimensional added wave resistance