the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
VISIR2: ship weather routing in Python
Gianandrea Mannarini
Mario Leonardo Salinas
Lorenzo Carelli
Nicola Petacco
Josip Orović
Ship weather routing, which involves suggesting lowemission routes, holds potential for contributing to the decarbonisation of maritime transport. However, including because of a lack of readily deployable opensource and openlanguage computational models, its quantitative impact has been explored only to a limited extent.
As a response, the graphsearch VISIR (discoVerIng Safe and effIcient Routes) model has been refactored in Python, incorporating novel features. For motor vessels, the angle of attack of waves has been considered, while for sailboats the combined effects of wind and sea currents are now accounted for. The velocity composition with currents has been refined, now encompassing leeway as well. Provided that the performance curve is available, no restrictions are imposed on the vessel type. A cartographic projection has been introduced. The graph edges are quickly screened for coast intersection via a Kdimensional tree. A leastCO_{2} algorithm in the presence of dynamic graph edge weights has been implemented and validated, proving a quasilinear computational performance. The software suite's modularity has been significantly improved, alongside a thorough validation against various benchmarks. For the visualisation of the dynamic environmental fields along the route, isochronebounded sectors have been introduced.
The resulting VISIR2 model has been employed in numerical experiments within the Mediterranean Sea for the entirety of 2022, utilising meteooceanographic analysis fields. For a 125 m long ferry, the percentage saving of overall CO_{2} expenditure follows a biexponential distribution. Routes with a carbon dioxide saving of at least 2 % with respect to the leastdistance route were found for prevailing beam or head seas. Twodigit savings, up to 49 %, were possible for about 10 d in a year. In the case of an 11 m sailboat, time savings increased with the extent of path elongation, particularly during upwind sailing. The sailboat's routes were made approximately 2.4 % faster due to optimisation, with the potential for an additional 0.8 % in savings by factoring in currents.
VISIR2 serves as an integrative model, uniting expertise from meteorology, oceanography, ocean engineering, and computer science, to evaluate the influence of ship routing on decarbonisation efforts within the shipping industry.
 Article
(8420 KB)  Fulltext XML
 Companion paper 1
 Companion paper 2

Supplement
(12136 KB)  BibTeX
 EndNote
As climate change, with its unambiguous attribution to anthropogenic activities, rapidly unfolds (IPCC, 2023), the causal roles played by various sectors of the economy, as well as the possibilities for mitigation, are becoming more evident. This holds true for the shipping sector as well (IPCC, 2022), which has begun taking steps to reduce its carbon footprint. The International Maritime Organization (IMO) adopted an initial decarbonisation strategy in 2018 (IMO, 2018a), which was later revised in 2023. The new ambition is to achieve complete decarbonisation by midcentury, addressing all greenhouse gas (GHG) emissions, with a partial uptake of nearzero GHG technology as early as 2030 (IMO, 2023). While no new measures have been adopted yet, the revised strategy is expected to boost efforts to increase the energy efficiency of shipping in the short term (Smith and Shaw, 2023). In line with the European Green Deal, the European Union has adopted new rules to include various GHG (carbon dioxide, methane, and nitrous oxide) emissions from shipping in its Emissions Trading System (EUETS), starting from 2024 (https://climate.ec.europa.eu/euaction/transportemissions/reducingemissionsshippingsector_en, last access: 20 May 2024). For the first time ever, this will entail surrendering allowances for GHG emissions from vessels as well.
Zeroemission bunker fuels are projected to cost significantly more than presentday fossil fuels (AlAboosi et al., 2021; Svanberg et al., 2018). Thus, minimising their use will be crucial for financial sustainability. This necessitates energy savings through efficient use, achieved via both technical (e.g. windassisted propulsion systems or WAPSs) and operational measures (e.g. speed reduction and ship weather routing). According to the “CEShip” model, a GHG emissions model for the shipping sector, a global reduction in GHG emissions by 2030 by up to 47 % relative to 2008 levels could be feasible through a combination of operational measures, technical innovations, and the use of nearzeroGHG fuels (Faber et al., 2023). A separate study focusing on the European fleet estimates that a reduction in sailing speed alone could potentially lead to a 4 %–27 % emission reduction, while combining technical and operational measures might provide an additional 3 %–28 % reduction (Bullock et al., 2020). The impact of speed optimisation on emissions varies significantly, with potential percentage savings ranging from a median of 20 % to as high as 80 %, depending on factors such as the actual route, meteorological and marine conditions, and the vessel type (Bouman et al., 2017). In that paper, while cases are reported of up to 50 % savings, the role of weather routing was generally assessed to be lower than 10 %.
The variability in percentage savings reported in the literature can be attributed to the diversity of routes considered, the specific weather conditions, and the types of vessels analysed. Additionally, reviews often use a wide range of bibliographic sources, including grey literature, company technical reports, white papers, and works that fail to address the actual meteorological and marine conditions.
The VISIR (discoVerIng Safe and effIcient Routes, pronunciation: /vi′zi:r/) ship weather routing model was designed to objectively assess the potential impact of oceanographic and meteorological information on the safety and efficiency of navigation. So far, two versions of the model have been released (VISIR1.a – Mannarini et al., 2016a, and VISIR1.b – Mannarini and Carelli, 2019a). However, the use of a proprietary coding language (MATLAB) may hinder its further adoption. Also, the experience with VISIR1 suggested the need to enhance the modularity of the package and implement other best practices of scientific coding, as recommended in Wilson et al. (2014). Another area where innovation seemed possible was the development of a comprehensive framework to perform weather routing for both motor vessels and sailboats. While some aspects of this requirement were covered through a more modular approach, it also involved rethinking how to utilise environmental fields such as waves, currents, and wind. Furthermore, while the carbon intensity savings of leasttime routes were already estimated via VISIR1.b, a dedicated algorithm for the direct computation of leastCO_{2} routes was lacking.
To address all these requirements, we designed, coded, tested, and conducted extensive numerical experiments with the VISIR2 model (Salinas et al., 2024a). VISIR2 is a Pythoncoded software, inheriting from VISIR1 the fact that it is based on a graphsearch method. However, VISIR2 is a completely new model, leveraging the previous experience while also offering many new solutions and capabilities. Part of the validation process made use of its ancestor model VISIR1. Computational performance has been enhanced, and efforts have been made to improve usability. In view of use for routing of WAPS vessels, we have developed a unified modelling framework, accounting for involuntary speed loss in waves, advection via ocean currents, thrust, and leeway from wind. A novel algorithm for finding the leastCO_{2} routes has been devised, building upon Dijkstra's original method. The visualisation of dynamic information along the route has been further improved using isochronebounded sectors. VISIR2 features are thoroughly described in this paper, along with some case studies and suggestions for possible development lines in the future.
The VISIR2 model enabled systematic assessment of CO_{2} savings deriving from optimal use of meteooceanographic information in the voyage planning process. Differently from most of the existing works cited above, the emissions could clearly be related to the significant wave height, relative angle of waves, and also engine load factor. This is discussed further in Sect. 6.1. Ocean currents can be added to the optimisation, providing more realistic figures for the emission savings.
VISIR2 is a numerical model wellsuited for both academic and technical communities and facilitating the exploration and quantification of ship routing's potential for operational decarbonisation. In the European Union, there may be a need for independent verification of emissions reported in the monitoring, reporting, and verifying (MRV) system (https://www.emsa.europa.eu/thetismrv.html, last access: 20 May 2024). Utilising baseline emissions computed through a weather routing model could enhance the accuracy and reliability of this verification process. Additionally, the incorporation of shipping into the EUETS intensifies the importance of minimising GHG emissions. Internationally, even with the adoption of zero or lowcarbon fuels encouraged by the recent update of the IMO's strategy, it remains critical to save these fuels as much as possible. Therefore, VISIR2 could potentially serve as a valuable tool for saving both fuel and allowances.
The remainder of this paper includes a literature investigation in Sect. 1.1 followed by an indepth presentation of the innovations introduced by VISIR2 in Sect. 2. Subsequently, the model validation is discussed in Sect. 3 and its performance is assessed in Sect. 4. Several case studies in the Mediterranean Sea follow (Sect. 5). The results of the CO_{2} emission percentage savings, some potential use cases, and an outlook on possible developments of VISIR2 are discussed in Sect. 6. Finally, the conclusions are presented in Sect. 7. The Appendix contains technical information regarding the computation of the angle of attack between ships' heading and course (Appendix A), as well as details about the neural network employed to identify the vessel performance curves (Appendix B).
1.1 Literature review on weather routing
This compact review of systems for ship weather routing will be limited to web applications (Sect. 1.1.1) and peerreviewed research papers (Sect. 1.1.2). It is further restricted to the freely available versions of desktop applications, while the selection of papers is meant to update the wider reviews already provided in Mannarini et al. (2016a) and Mannarini and Carelli (2019b). A critical gap analysis (Sect. 1.1.3) completes this subsection.
1.1.1 Web applications
FastSeas (https://fastseas.com/, last access: 20 May 2024) is a weather routing tool for sailboats, with editable polars and the possibility of considering a motor propulsion. Wind forecasts are taken from the National Oceanic and Atmospheric Administration Global Forecast System (NOAA GFS) model and ocean surface currents from NOAA Ocean Surface Current Analyses Realtime (OSCAR). The tool makes use of Windy imagery, a free choice of endpoints is offered, departure time can vary between the present and a few days in the future, and the voyage plan is exportable in various formats.
The Avalon web router (https://www.webrouter.avalonrouting.com/computeroute, last access: 20 May 2024) provides a coastal weather routing service for sailboats within subregions of France, the United Kingdom, and the United States. It offers a choice among tens of sailboats, with the option to also consider a motorassisted propulsion. Hourly departure dates within a couple of days are allowed, and ocean weather along the routes is provided in tabular form.
GUTTAVISIR (https://guttavisir.eu, last access: 20 May 2024) is a weather routing tool for a specific ferry developed in the frame of the “savinG fUel and emissions from mariTime Transport in the Adriatic region” (GUTTA) project. It provides both leasttime and leastCO_{2} routes between several ports of call. It makes use of operational wave and current forecast fields from the Copernicus Marine Environment Monitoring Service (CMEMS). The route departure date and time and the engine load can be varied. Waves, currents, or isolines can be rendered along with the routes, which can be exported.
OpenCPN (https://opencpn.org/, last access: 20 May 2024) is a comprehensive opensource platform, including a weather routing tool for sailboats. The product name originates from “Opensource Chart Plotter Navigation”. The vessel performance curves are represented via polars, and forecast data in GRIB format or a climatology can be used. Nautical charts can be downloaded and integrated into the graphical user interface. The programming language is C++, and a velocity composition with currents is accounted for. A detailed documentation of the numerical methods used is lacking though.
1.1.2 Research papers
A review of ship weather routing methods and applications was provided by Zis et al. (2020). Several routing methods such as the isochrone method, dynamic programming, calculus of variations, and pathfinding algorithms were summarised before a taxonomy of related literature was proposed. The authors made the point that the wide range of emission savings reported in the literature might in future be constrained via defining a baseline case; providing benchmark instances; and performing sensitivity analyses, e.g. on the resolution of the environmental data used.
A specific weather routing model was documented by Vettor and Guedes Soares (2016). It integrates advanced seakeeping and propulsion modelling with multiobjective (fuel consumption, duration, and safety) optimisation. An evolutionary algorithm was used, with initialisation from both random solutions and singleobjective routes from a modified Dijkstra's algorithm. Safety constraints were considered via probabilities of exceeding thresholds for slamming, green water, or vertical acceleration. A specific strategy was proposed to rank solutions within a Pareto frontier.
A stochastic routing problem was addressed in Tagliaferri et al. (2014). A single upwind leg of a yacht race is considered, with the wind direction being the stochastic variable. The vessel was represented in terms of polars, and the optimal route was computed via dynamic programming. The skipper's risk propensity was modelled via a specific preference on the wind transition matrices. This way it was shown how the risk attitude affects the chance to win a race.
Ladany and Levi (2017) developed a dynamicprogramming approach to sailboat routing which accounts for the tacking time. The latter was assumed to be proportional to the amplitude of the course change. Furthermore, a sensitivity analysis was conducted, considering both the uncertainty in wind direction and the magnitude of the discretisation step in the numerical solution.
Sidoti et al. (2023) provided a consistent framework for a dynamicprogramming approach for sailboats, considering both leeway and currents. In order to constrain the course of the boat on the available edges on the numerical grid, an iterative scheme was adopted. Case studies with NAVGEM winds and global HYCOM currents were carried out in a region across the Gulf Stream. The results without leeway were validated versus OpenCPN.
The impact of stochastic uncertainty on WAPS ships was addressed by Mason et al. (2023b). A dynamicprogramming technique was used, and “a priori” routing (whereby information available at the start of the journey only is used) was distinguished from “adaptive” routing (whereby the optimum solution is updated based on information that becomes available every 24 h along a journey). The latter strategy is shown, for voyages lasting several days, to be more robust with respect to the unavoidable stochastic uncertainty in the forecasts.
1.1.3 Knowledge gap
A few open web applications exist, mainly for sailboats and with limited insight into their numerical methods. Case study results from weather routing systems developed in academia have been published, but (with the exception of Mason et al., 2023b) no systematic assessment of CO_{2} savings has been provided. Furthermore, no related software has been disclosed in any case. The OpenCPN package lacks both a peerreviewed publication and documentation of the methods implemented. A prevalence of dynamicprogramming approaches is noted, especially for web applications, with the graphsearch method being used in research papers only. The tools focus either on sailboats (with or without a motor) or on motor vessels. When both are available (such as in FastSeas), the motor vessel is described in terms of polars.
From this assessment, the lack of an opensource and welldocumented ship weather routing model, for both motor vessels and sailboats, with flexible characterisation of vessel performance appears as a gap which the present work aims to close.
This section includes a revision of the vessel kinematics of VISIR, as given in Sect. 2.1; changes in the graph generation procedure and in the use of static environmental fields in Sect. 2.2; updates to the computation of graph edge weights in Sect. 2.3; an additional optimisation objective in the shortestpath algorithm in Sect. 2.4; new vessel performance models in Sect. 2.5; innovative visualisation capabilities in Sect. 2.6; and a more modular structure of the software package, presented in Sect. 2.7. Further technical details of the VISIR2 code are presented in the software manual, provided along with its online repository (Salinas et al., 2024a).
2.1 Kinematics
For a graphsearch model such as VISIR to deal with waves, currents, and wind for various vessel types, several updates to its approach for velocity composition were needed. They included both generalisations and use of new quantities, addressed in this subsection, as well as a new numerical solution, addressed in Appendix A.
As in Mannarini and Carelli (2019b), the kinematics of VISIR2 is based on both the principle of superposition of velocities and the discretised sailing directions existing on a graph. However, while previously the vessel's speed over ground (SOG) was obtained from the magnitude of the vector sum of speed through water (T = STW $\widehat{t}$) and ocean current (w), we show here that, more generally, SOG is the magnitude of a vector G = SOG$\phantom{\rule{0.125em}{0ex}}\widehat{e}$, being the sum of the forward speed F = $F\widehat{h}$ and an effective current ω, and both will be defined in the following. In the absence of leeway, the F and T vectors are identical and ω=w, so the newer approach encompasses the previous one.
Making reference to Fig. 1, the use of a graph constrains the vessel's course over ground to be along $\widehat{e}$, being the orientation of one of the graph's arcs. Thus, any cross component of velocity or along a $\widehat{o}$ versor such that $\widehat{e}\cdot \widehat{o}=\mathrm{0}$ must be null.
This implies that, to balance the cross flow from the currents, the vessel must head into a direction $\widehat{h}$ slightly different from the course $\widehat{e}$. In Mannarini and Carelli (2019b), such an angle of attack was defined as
where for both ψ_{h} (heading, or HDG) and ψ_{e} (course over ground, or COG) a nautical convention is used (due north, to direction). That framework is here generalised to also deal with the vector nature of some environmental fields, such as waves or wind. Using a meteorological convention (due north, from direction) for both the ψ_{a} (waves) and ψ_{i} (wind) directions, we here introduce the δ_{f} angles, defined as
where ${\mathit{\gamma}}_{f}={\mathit{\psi}}_{f}{\mathit{\psi}}_{e}$ constitutes the relative angle between the f environmental field and the ship's course. f=a for waves and f=i for wind. In computing angular differences, their convention should be considered (see the angular_utils.py
function in the VISIR2 code). Thus, δ_{f}=0 whenever the ship heads in the direction from which the field comes, and γ_{f}=0 if her course is in such a direction.
Furthermore, we here define leeway as a motion, caused by the wind, transversal to the ship's heading. From the geometry shown in Fig. 1, the oriented direction of leeway ψ_{L} is given by
with
where the ⌊⋅⌋ delimiters indicate the floor function. Thus, ϵ is positive for δ_{i} in the [0,180°] range and flips every 180° of the argument. This is not the sole possible definition of leeway. For instance, in Breivik and Allen (2008) a distinction between the downwind and crosswind component of leeway is made. However, the present definition is consistent with the subsequent data of vessel performance (Sect. 2.5.2).
Upon expressing the magnitudes F of the vessel's forward velocity and L of leeway velocity as
with wind magnitude V_{i}, significant wave height H_{s}, and engine load χ, a leeway angle, being the ship's heading change between the F and T vectors, can be defined as
The aboveintroduced α_{L} is not a constant but depends on both wind magnitude V_{i} and the module of the relative direction, $\left{\mathit{\delta}}_{i}\right$. Also, from Fig. 1, it is seen that F = STW ⋅cos α_{L}. Thus, in the absence of leeway (α_{L}=0), one retrieves the identity of F and STW, which was an implicit assumption made in Mannarini and Carelli (2019b).
Equations (5a) and (5b) comprise a major innovation with respect to the formalism of Mannarini and Carelli (2019b), as an angular dependence in the performance curve is also introduced in VISIR for motor vessels for the first time (Sect. 2.5). Equation (5a) includes dependencies on both wind and waves. Furthermore, a possible dependency on χ, the fractional engine load (or any other propulsion parameter), is highlighted here.
Within this formalism, if the vessel is a sailboat (or rather a motor vessel making use of a WAPS), just one additional condition should be considered. That is, given the windmagnitudedependent nogo angle α_{0}(V_{i}), upwind navigation is not permitted, or
Now, given a water flow expressed by the vector
and making reference to Fig. 1, the flow projections along ($\widehat{e}$) and across ($\widehat{o}$) the vehicle course are, respectively,
where for the ocean flow direction ψ_{w} the nautical convention is also used.
In analogy to Eqs. (9a) and (9b) and also using nautical convention for ψ_{L}, the along and crosscourse projections of the leeway are given by
The simple relations on the righthand side (r.h.s.) of Eqs. (10a) and (10b) follow from the similitude of the red and greenshaded triangles in Fig. 1. As δ is typically a small angle (cf. Appendix A), it is apparent that the cross component of the leeway, ${w}_{\perp}^{\left(L\right)}$, is the dominant one. Its sign is such that it is always downwind; see Fig. 1. If relevant, the Stokes drift (van den Bremer and Breivik, 2018) could be treated akin to an ocean current, and one would obtain for its projections a couple of equations formally identical to Eqs. (9a) and (9b).
Finally, the components of the effective flow ω advecting the vessel are
Due to Eqs. (10a) and (10b), both ω_{∥} and ω_{⊥} are functions of δ. We recall that the “cross” and “along” specifications refer to vessel course ψ_{e}, differing from vessel heading by the δ angle.
The graphical construction in Fig. 1 makes it clear that G = SOG $\widehat{e}$ equals T+w or F+ω. Using the latter equality, together with the course assignment condition, and projecting along both $\widehat{e}$ and $\widehat{o}$, two scalar equations are obtained, namely
Equations (12a) and (12b) are formally identical to those found in Mannarini and Carelli (2019b) in the presence of ocean currents only. This fact supports the interpretation of ω as an effective current.
However, Eqs. (12a) and (12b) alone are no more sufficient to determine the ocean current vector w. In fact, it is mingled with the effect of wind through leeway to form the effective flow ω (Eqs. 11a, 11b). This is why, in the presence of strong winds, reconstruction of ocean currents from data of COG and HDG via a naive inversion of Eqs. (12a) and (12b) is challenging. This was indeed found by Le Goff et al. (2021) using automatic identification system (AIS) data across the Agulhas Current.
As it reduces the ship's speed available along its course (Eq. 12a), the angle of attack δ plays a pivotal role in determining SOG. However, in the presence of an angledependent vessel speed (Eq. 5a), δ is no longer given by a simple algebraic equation corresponding to Eq. (12b) as in Mannarini and Carelli (2019b) but by a transcendental one:
In fact, due to Eq. (2), the r.h.s. of Eq. (13) depends both explicitly and implicitly on δ. Just in the limiting case of null currents, Eqs. (6), (10b), and (12b) collectively imply that
However, in general, the actual value of the forward speed F is only determined once the δ angle is retrieved from Eq. (13). To our knowledge, the transcendental nature of the equation defining δ had not been pointed out previously. Furthermore, in VISIR2 an efficient numerical approximation for its solution is provided; see Appendix A. This is also a novelty with practical benefits for the development of efficient operational services based on VISIR2.
We note that Eq. (13) holds if and only if
Should this not be the case, the vessels' forward speed would not balance the effective drift.
As F is always nonnegative, Eq. (13) implies that $\mathrm{sgn}\left(\mathit{\delta}\right)=\mathrm{sgn}\left({\mathit{\omega}}_{\perp}\right)$. In particular, in the case of an effective cross flow ω_{⊥} bearing, as in Fig. 1, to starboard, an anticlockwise change in vehicle heading (δ<0) is needed for keeping course (note $\widehat{o}$ versor's orientation in Fig. 1).
Equations (12a) and (12b) can be solved for the speed over ground (SOG), which reads
According to Eq. (16), the cross flow ω_{⊥} always reduces SOG, as part of vehicle momentum must be spent balancing the drift. The alongedge flow ω_{∥} (or “effective drag”) may instead either increase or decrease SOG.
Finally, given that $\mathit{G}=\mathrm{d}\mathit{x}/\mathrm{d}t$, by taking the module of the left side, and approximating the r.h.s. with its finitedifference quotient, the graph edge weight δt is computed as
where δx is the edge length and SOG is given by Eq. (16). As the environmental fields determining SOG are both space and time dependent, the weights δt are computed via the specific interpolation procedures in Sect. 2.3 and the shortest paths via the algorithms provided in Sect. 2.4.
From Eq. (17), it follows that the condition
should be checked in case the specific graphsearch method used does not allow for use of negative edge weights (as is the case for Dijkstra's algorithm; cf. Bertsekas, 1998). Violation of Eq. (18) may occur in the presence of a strong counterflow ω_{∥} along a specific graph edge, whose weight must correspondingly be set to not a number.
The CO_{2} emissions along a graph edge are given by
where $\mathrm{\Gamma}=\mathrm{\Gamma}\left(\right{\mathit{\delta}}_{i},{V}_{i};{\mathit{\delta}}_{a},{H}_{\mathrm{s}},\mathit{\chi})$ is the CO_{2} emission rate of the specific vessel in the presence of the actual meteomarine conditions. Both δt and δCO_{2} are used in the leastCO_{2} algorithm introduced in Sect. 2.4.
2.2 Graph generation
Graph preparation is crucial for any graphsearch method. Indeed, graph edges represent potential route legs. Therefore, the specific edges included within the graph directly influence the route topology. In addition, as will be shown in Sect. 2.3.2, the length of edges affects the environmental field value represented by each edge. On the other hand, graph nodes determine the accessible locations within the domain. Therefore, they must be selected with consideration of both the presence of landmasses and shallow waters.
The structure of the mesh is also the most fundamental difference between a graphsearch method (such as Dijkstra's or A^{*}) and dynamic programming. Indeed, a dynamicprogramming problem can be transformed into a shortestpath problem on a graph (Bertsekas, 1998, Sect. 2.1). In the former, the nodes are organised along sections of variable amplitude corresponding to stages. In the latter, nodes can uniformly cover the entire domain. However, for both dynamic programming (Mason et al., 2023b) and graph methods (Mannarini et al., 2019b), the number of edges significantly impacts the computational cost of computing a shortest path.
To efficiently address all these aspects, VISIR2's graph preparation incorporates several updates compared to its predecessor, VISIR1b, as outlined in the following subsections.
2.2.1 Cartographic projection
The nautical navigation purpose of a ship routing model necessitates the provision of both length and course information for each route leg. On the other hand, environmental fields used to compute a ship's SOG (Eq. 16) are typically provided as a matrix of values in spherical coordinates (latitude and longitude, similar to the fields used in Sect. 5.1). To address these aspects, VISIR2 introduces a projection feature, which was overlooked in both VISIR1.a and VISIR1.b. The Mercator projection was chosen for its capacity to depict constant bearing lines as straight lines (Feeman, 2002). The implementation was carried out using the pyproj
library, which was employed to convert spherical coordinates on the WGS 84 (World Geodetic System 1984) ellipsoid into a Mercator projection, with the Equator serving as the standard parallel. For visualisation purposes, both the Matplotlib
and cartopy
libraries were utilised.
2.2.2 Edge direction
Leg courses of a ship route originate from graph edge directions. To determine them, we consider the Cartesian components of the edges in a projected space. As seen from Fig. 2, this approach results in smaller angles relative to due north compared to using an unprojected graph. For the graph (ν, Δx) values and average latitude of the case studies considered in this paper (Sect. 5), the maximum error between an unprojected graph and a Mercator projection would be about 5° (cf. Table 1). However, this angle critically impacts vessels like sailboats, whose performance hinges on environmental field orientation. Additionally, at higher latitudes, the courses tend to cluster along meridional directions. To achieve a more isotropic representation of courses, constructing the graph as a regular mesh in the projected space would be needed, which is left for future refinements of VISIR2.
For a given sea domain, a graph is typically computed once and subsequently utilised for numerous different routes. However, in VISIR1, the edge direction was recalculated every time a graph was utilised. In VISIR2, the edge direction is computed just once, during graph generation, after which it is added to a companion file of the graph (edge_orientation
). In the definition of the edge direction, the nautical convention (due north, to direction) is used in VISIR2 graphs.
Finally it should be noted that in a directed graph, such as VISIR2's, edge orientation also matters. Each edge acts as an arrow, conveying flow from “tail” to “head” nodes. Edge orientation refers to the assignment of the head or tail of an edge. Orientation holds a key to understanding vessel performance curves, explored further in Sect. 2.5.
2.2.3 Quasicollinear edges
A further innovation regarding the graph involves an edge pruning procedure. This eliminates redundant edges that carry very similar information. Indeed, some edges have the same ratio of horizontal to vertical grid hops (see the possible_hops()
function). Examples of these edges are the solid and dashed arrows in Fig. 2b. When the grid step size and the number of hops are not too large, these edges point to nearly the same angle relative to north. However, starting from an equidistant grid in spherical coordinates, the vertical spacing in Mercator projection is uneven. Thus, the directions of those edges are not exactly the same, and we call such edges “quasicollinear”. In VISIR2, only the shortest one among these quasicollinear edges is retained. This corresponds to the solid arrows in Fig. 2b. This reduces the number N_{q1} of edges within a single quadrant to
where φ is Euler's totient function and ν the maximum number of hops from a given node of the graph. Thus, the quantity right of the inequality represents the total number of edges of a quadrant, including quasicollinear ones. Using Eq. (20), at ν=4 already more than onethird of all edges are pruned and, at ν=10, nearly half of them are (cf. Table 1).
This benefits both the computer memory allocation and the computing time for the shortest path. The latter is linear in the number of edges; cf. Bertsekas (1998, Sect. 2.4.5). A further benefit of pruning quasicollinear edges is a more faithful representation of the environmental conditions. In fact, the environmental field's values at the graph nodes are used for estimating the edge weights (see Sect. 2.3). Thus, keeping just shorter edges avoids using less spatially resolved information.
2.2.4 Bathymetry and draught
The minimal safety requirement is that navigation does not occur in shallow waters. This corresponds to the condition that vessel draught T does not exceed sea depth z at all graph edges used for the route computations. This is equivalent to a positive under keel clearance (UKC) of z−T. As explained later in Sect. 2.3.2, this can be checked by either evaluating the average UKC at the two edge nodes or interpolating it at the edge barycentre.
However, for a specific edge, UKC could still be positive and the edge cross the shoreline. This is avoided in VISIR by checking for mutual edge–shoreline crossings. Given the burden of this process, in VISIR1b a procedure for restricting the check to inshore edges was introduced. In VISIR2, as envisioned in Mannarini and Carelli (2019b, Appendix C), the process of searching for intersections is carried out using a Kdimensional tree (KDT; Bentley, 1975; Maneewongvatana and Mount, 1999). This is a means of indexing the graph edges via a spatial data structure which can effectively be queried for both nearest neighbours (coast proximity of nodes) and range queries (coast intersection of edges). The scipy.spatial.KDTree
implementation was used (https://docs.scipy.org/doc/scipy/reference/, last access: 20 May 2024). Use of an advanced data structure in the graph is also a novelty of VISIR2.
Various bathymetric databases can be used by VISIR2. For European seas, the EMODnet dataset (https://portal.emodnetbathymetry.eu/, last access: 20 May 2024) ($\mathrm{1}/\mathrm{16}$ arcmin resolution or about 116 m) was used while, for global coverage, GEBCO_2022 (https://www.gebco.net/data_and_products/gridded_bathymetry_data/, last access: 20 May 2024) (15 arcsec resolution or about 463 m) is available.
2.2.5 Shoreline
The bathymetry dataset, if detailed enough, can even be used for deriving an approximation of the shoreline. From Fig. 3 it is seen that a “pseudoshoreline” derived from the UKC = 0 contour line of a bathymetry that is fine enough (the EMODnet one) can effectively approximate an official shoreline (the GSHHG (https://www.ngdc.noaa.gov/mgg/shorelines/, last access: 20 May 2024) one, at the “high”type resolution of 200 m).
Such a pseudoshoreline is the one used in VISIR2 for checking the edge crossing condition specified in Sect. 2.2.4.
2.3 Edge weights
For computing shortest paths on a graph, its edge weights are preliminarily needed. Due to Eqs. (16) and (5a)–(5b), they depend on both space and timedependent environmental fields, whose information has to be remapped to the numerical grids of VISIR2. This is done in a partly differently way than in VISIR1, providing users with improved flexibility and control over numerical fields. These novel options are documented in what follows.
2.3.1 Temporal interpolation
The time at which edge weights are evaluated is key to the outcome of the routing algorithm. In Mannarini et al. (2019b), to improve on the coarse time resolution of the environmental field, a linear interpolation in time of the edge weight was introduced (“Tint =1” option in Fig. 4). In VISIR2, instead, the environmental field values (grey dots) are preliminarily interpolated in time on a finer grid with Δτ spacing (“Tint =2” or blue dots). Then, the edge weight at the nearest available time step (np.floor
function used, corresponding to the blue segments) is selected.
2.3.2 Spatial interpolation
The numerical environmental field and the graph grid may possess varying resolutions and projections. Even if they were identical, the grid nodes might still be staggered. Moreover, it is necessary to establish a method for assigning the field values to the graph edges. For all these reasons, spatial interpolation of the environmental fields is essential.
VISIR2 first computes an interpolant using the scipy.interpolate.interp2d
method. Then, to assign a representative value of the φ field on an edge, two options are available, which are depicted in Fig. 5. In the first one or “Sint =0”, an average between the edge head and tail values is computed, ${\mathit{\phi}}_{\mathrm{0}}=({\mathit{\phi}}_{h}+{\mathit{\phi}}_{t})/\mathrm{2}$. In the second one or “Sint =1”, the field interpolator is evaluated at the location of the edge barycentre, ${\mathit{\phi}}_{\mathrm{1}}=\mathit{\phi}\left(({\mathit{x}}_{h}+{\mathit{x}}_{t}\right)/\mathrm{2})$. As the field is generally nonlinear, this leads to different outcomes.
The two interpolation schemes were tested with various nonconvex functions, and some results are reported in Sect. S0 of the Supplement. We observed that both options converge to a common value as the $(\mathrm{1}/\mathrm{\Delta}x)$ resolution of the graph grid increases. Thus, both options benefit from pruning collinear edges (Sect. 2.2.3), as they remove some longer edges from the graph, but neither demonstrates consistent superiority over the other in terms of fidelity.
However, setting Sint =0 results in higher computational efficiency. This is because the interpolator is applied at each node with Sint =0, whereas it is applied at each edge with Sint =1. Given that the number of edges exceeds the number of nodes by a factor defined by Eq. (20), for the case studies ($\mathrm{4}\le \mathit{\nu}\le \mathrm{5}$), the computational time for Sint =0 was approximately 1 order of magnitude smaller than with Sint =1. Therefore, the latter was chosen as the default for VISIR2.
The resulting edgerepresentative environmental field value affects the edge delay (Eq. 17), consequently impacting the vessel's speed as per Eqs. (5a) and (5b). Hence, a nonlinear relationship exists between the field value and the local sailing speed. Thus, the actual choice of the Sint parameter is not expected to systematically bias the vessel's speed in either direction.
Even before the spatial interpolation is performed, the socalled “seaoverland” extrapolation is applied to the marine fields. This step, which is needed for filling the gaps in the vicinity of the shoreline, is conceptually performed as in Mannarini et al. (2016a, Fig. 7) and implemented in the seaoverland
function.
Since wave direction is a circular periodic quantity, we calculate its average using the circular mean (https://en.wikipedia.org/wiki/Circular_mean, last access: 20 May 2024) ahead of interpolation. This differs from wind direction, typically given as Cartesian components (u,v), which can be interpolated directly.
2.4 Shortestpath algorithms
A major improvement made possible by the Python coding of VISIR2 is the availability of builtin, advanced data structures such as dictionaries, queues, and heaps. They are key in the efficient implementations of graphsearch algorithms (Bertsekas, 1998). In particular, as data structures are used, Dijkstra's algorithm worstcase performance can improve from quadratic, 𝒪(N^{2}), to linear–logarithmic, 𝒪(Nlog N), where N is the number of graph nodes.
Nonetheless, Dijkstra's original algorithm exclusively accounted for static edge weights (Dijkstra, 1959). When dynamic edge weights are present, Orda and Rom (1990) demonstrated that, in general, there are no computationally efficient algorithms. However, they also showed that, upon incorporating a waiting time at the source node, it is possible to keep the algorithmic complexity of a static problem. Given the assumption that the rate of variation in the edge delay δt satisfies
an initial waiting time is not even unnecessary. This condition was assumed to hold in Mannarini et al. (2016a) for implementing a timedependent Dijkstra's algorithm. That version of the shortestpath algorithm could not be used with an optimisation objective differing from voyage duration. As one aims to compute, for example, leastCO_{2} routes, the algorithm requires further generalisation. This has been addressed in VISIR2 via the pseudocode provided in both Algorithms 1 and 2. For its implementation in Python, we made use of a modified version of the single_source_Dijkstra
function of the NetworkX
Python library. The modification consisted in retrieving an edge weight at a specific time step. This is achieved via Algorithm 2. In this, the cost.at_time
pseudofunction represents a NetworkX
method to access the edge weight information.
We note the generality of the pseudocode with respect to the edge weight type (wT parameter in both Algorithms 1 and 2). This implies that the same algorithm could also be used to compute routes that minimise figures of merit differing from CO_{2}, such as different GHG emissions, cumulated passenger comfort indexes, or total amount of underwater radiated noise. This hinges solely on the availability of sufficient physical–chemical information to compute the corresponding edge weight in relation to the environmental conditions experienced by the vessel. A flexible optimisation algorithm is yet another novelty of VISIR2.
The shortestdistance and the leasttime algorithms invoked for both motor vessels and sailboats are identical. Differences occur at the postprocessing level only, as different dynamical quantities (related to the marine conditions or the vessel kinematics) have to be evaluated along the optimal paths. Corresponding performance differences are evaluated in Sect. S1 of the Supplement.
NetworkX
graph, source and target nodes, type of edge weight, maximum number of time steps, and time resolution, respectively2.4.1 NonFIFO
As previously mentioned, Dijkstra's algorithm can recover the optimal path in the presence of dynamic edge weights if Eq. (21) is satisfied. For cases where the condition is not met, Orda and Rom (1990) presented a constructive method for determining a waiting time at the source node. In this case, waiting involves encountering more favourable edge delay values, leading to the computation of a faster path. In other words, employing a firstin–firstout (FIFO) strategy for traversing graph edges may not always be optimal. This is why such a scenario is referred to as “nonFIFO”. We observe that condition Eq. (21) is violated when a graph edge, which was initially unavailable for navigation, suddenly becomes accessible. While this is a rather infrequent event for motor vessels (in Mannarini and Carelli, 2019b, nonFIFO edges were just 10^{−6} of all graph edges), it is a more common situation for sailboats navigating in areas where the wind is initially too weak or within the nogo zone (Eq. 7, Fig. 7). Indeed, the unavailability of an edge can be suddenly lifted as the wind strengthens or changes direction. However, under a FIFO hypothesis, the leasttime algorithm would not wait for this improvement of the edge delay to occur. Rather, it would look for an alternative path that avoids the forbidden edge, potentially leading to a suboptimal path. In the case study of this paper, such a situation occurred for about $\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$ of the total number of sailboat routes; see Sects. 2.5.2 and S3.2.
2.5 Vessel modelling
At the heart of the VISIR2 kinematics of Sect. 2.1 are the vessel forward and transversal speed in a seaway, Eqs. (5a)–(5b). In what follows, such a vessel performance function is also termed a “vessel model”.
In VISIR1 the forward speed resulted, for motor vessels, from a semiempirical parameterisation of resistances (Mannarini et al., 2016a) and, for sailboats, from polar diagrams (Mannarini et al., 2015). The transversal speed due to leeway was neglected.
In VISIR2 new vessel models were used, and just two of them are presented in this paper: a ferry and a sailboat. However, any other vessel type can be considered, provided that the corresponding performance curve is utilised in the Navi
module (Table 4, Fig. 8). The computational methods used to represent both the ferry and the sailboats are briefly described in Sect. 2.5.1–2.5.2. All methods provide the relevant kinematic quantities and, where applicable, the emission rates in correspondence to discrete values of the environmental variables. Such a “lookup table” (LUT) was then interpolated to provide VISIR2 with a function to be evaluated in the actual environmental (wave, currents, or wind) conditions. Additional LUTs can be used as well, and the relevant function for this part of the processing is vesselModels_identify.py
.
The interpolating function was either a cubic spline (for sailboats) or the outcome of a neuralnetworkbased prediction scheme (for the ferry). The neural network features are provided in Appendix B. While the neural network generally demonstrated superior performance in fitting the LUTs (see Sect. S2 of the Supplement), it provided unreliable data in extrapolation mode, as shown in Figs. 6–7. In contrast, the spline, when extrapolation was requested, returned the value at the boundary of the input data range.
2.5.1 Ferry
The ferry modelled in VISIR2 was a mediumsize RoPax vessel whose parameters are reported in Table 2.
A vessel's seakeeping model was used at the ship simulator at the University of Zadar, as documented in Mannarini et al. (2021). Therein, additional details about both the simulator and the vessel can be found. The simulator applied a forcing from wind waves of significant wave height H_{s} related to the local wind intensity V_{i} by
This relationship was derived by Farkas et al. (2016) for the wave climate of the central Adriatic Sea. The simulator then recorded the resulting vessel speed, as well as some propulsion and emissionrelated quantities. Leeway could not be considered by the simulator. The postprocessed data feed the LUT to then be used for interpolating both STW and the CO_{2} emission rate Γ as functions of significant wave height H_{s}, relative windwave direction δ_{a}=δ_{i}, and fractional engine load χ. The results are displayed in Fig. 6.
In a given sea state, the sustained speed is determined by the parameter χ. For head seas (δ_{a}=0°) STW is seen to decrease with H_{s}. The maximum speed loss varies from about 45 % of the calm water speed at χ=1 to about 70 % at χ=0.7 (Fig. 6a). For χ=0.7, the STW sensitivity to H_{s} decreases from head (δ_{a}=0°) to following seas (δ_{a}=180°, Fig. 6b). For this specific vessel, the increase in roll motion in beam seas, as discussed in Guedes Soares (1990), and its subsequent impact on speed loss do not appear to constitute a relevant factor.
The Γ rate, which is on the order of 1 t CO_{2} h^{−1}, shows a shallow dependence on H_{s} (Fig. 6c), while it is much more critically influenced by both χ and δ_{a} (Fig. 6c, d).
2.5.2 Sailboat
Any sailboat described in terms of polars can in principle be used by VISIR2. For the sake of the case study, a Bénétau First 36.7 was considered. Its hull and rigging features are given in Table 3.
The modelling of the sailboat STW was carried out by means of the WinDesign velocity prediction program (VPP). The tool was documented in Claughton (1999, 2003) and references therein. The VPP is able to perform a fourdegreesoffreedom analysis, taking into account a wide range of semiempirical hydrodynamic and aerodynamic models. It solves an equilibrium problem by a modified multidimensional Newton–Raphson iteration scheme. The analysis considered the added resistance due to waves by means of the socalled “Delft method” based on the Delft Systematic Yacht Hull Series (DSYHS). Besides, the tool allows us to introduce response amplitude operators derived from other techniques as well, such as computational fluid dynamics. The windwave relationship was assumed to be given by Eq. (22). For each wind configuration (i.e. speed and direction), the optimal choice of sail set was considered. The main sail and the jib sail were considered for upwind conditions; otherwise the combination of main sail and spinnaker was used.
The outcome corresponds to Eqs. (5a) and (5b) and is provided in Fig. 7. The nogo angle α_{0} varies from 27 to 53° as the wind speed increases from 5 to 25 kn. At any true wind angle of attack δ_{i}, the forward speed F increases with wind intensity, especially at lower magnitudes (Fig. 7a). The peak boat speed is attained for broad reach (δ_{i}≈135°). Leeway magnitude L instead is at its largest for points of sail between the nogo zone (δ_{i}=α_{0}) and beam reach (δ_{i}=90°); see Fig. 7b. As the point of sail transitions from the nogo zone to running conditions, the leeway angle α_{L} gradually reduces from 6 to 0°. This decrease follows a roughly linear pattern, as depicted in Fig. 7c.
2.6 Visualisation
Further innovations brought in by VISIR2 concern the visualisation of the dynamic environmental fields and the use of isolines.
To provide dynamic information via a static picture, the fields are rendered via concentric shells originating at the departure location. The shape of these shells is defined by isochrones. These are lines joining all sea locations which can be reached from the origin upon sailing for a given amount of time. This way, the field is portrayed at the time step the vessel is supposed to be at that location. Isochrones bulge along gradients of a vessel's speed. Such shells represent an evolution of the stripewise rendering introduced in VISIR1.b (Mannarini and Carelli, 2019b, Fig. 5). The saved temporal dimensional of this plot type allows for its application in creating movies, where each frame corresponds to varying values of another variable, such as the departure date or engine load; see the “Video supplement” of this paper. This visual solution is another novelty introduced by VISIR2.
In addition to isochrones, lines of equal distance from the route's origin (or “isometres”) and lines of equal quantities of CO_{2} emissions (or “isopones”) are also computed. The name isopone is related to energy consumption (the Greek word means “equal effort”), which, for an internal combustion engine, the CO_{2} emission is proportional to. Isopones bulge against gradients of emissions. Isometres do not bulge, unless some obstruction (shoals, islands, landmass in general) prevents straight navigation. Given that rendering is on a Mercator map, “straight” refers to a ship's movement along a constant bearing line. Isochrones correspond to the reachability fronts used in a model based on the level set equation (LSE) by Lolla (2016).
2.7 Code modularity and portability
Software modularity has been greatly enhanced in VISIR2. While in VISIR1 modularity was limited to the graph preparation, which was detached from the main pipeline (Mannarini et al., 2016a, Fig. 8), the VISIR2 code is organised into many more software modules. Their characteristics are given in Table 4, and the overall model workflow is shown in Fig. 8.
The modules can be run independently and can optionally save their outputs. Through the immediate availability of products from previously executed modules, this favours research and development activities. For operational applications (such as GUTTAVISIR) instead, the computational workflow can be streamlined by avoiding the saving of the intermediate results. VISIR2 module names are Italian words. This is done for enhancing their distinctive capacity; cf. Wilson et al. (2014). More details on the individual modules can be found in the user manual, provided as part of the present release (Salinas et al., 2024a).
A preliminary graphical user interface (GUI) is also available. In the version of VISIR2 released here, it facilitates the ports' selection from the World Port Index (https://msi.nga.mil/Publications/WPI, last access: 20 May 2024) database.
VISIR2 was developed on macOS Ventura (13.x). However, both path parameterisation and the use of virtual environments ensure portability, which was successfully tested for both Ubuntu 22.04.1 LTS and Windows 11, on both personal computers and two distinct highperformance computing (HPC) facilities.
Validating a complex model like VISIR2 is imperative. The code was developed with specific runs of VISIR1 as a benchmark. The validation of VISIR1 involved comparing its outcomes with both analytical and numerical benchmarks and assessing its reliability through extensive utilisation in operational services (Mannarini et al., 2016b).
Previous studies have compared VISIR1 to analytical benchmarks for both static wave fields (“Cycloid”; Mannarini et al., 2016a) and dynamic currents (“Techy”; Mannarini and Carelli, 2019b). Here, we present the results of executing the same tests with VISIR2, at different graph resolutions, as shown in Table 5. The errors are consistently found to be below 1 %.
Additionally, in Mannarini et al. (2019b, Table II), routes computed in dynamic wave fields were compared to the results from a model based on the LSE. For these specific runs of VISIR2, the benchmarks were taken as LSE simulations at the nearest grid resolution. Notably, VISIR2 consistently produces shorterduration routes compared to VISIR1.b (see Table 6).
For both analytical and numerical benchmarks, distinct from the scenario discussed in Sect. 2.2.3, quasicollinear edges were retained in the graphs.
During the tests mentioned earlier, the vessel's STW remained unaffected by vector fields. In instances where there was a presence of current vector fields (Techy oracle), they were merely added to STW, without directly impacting it. Therefore, the enhanced capability of VISIR2 to accommodate angledependent vessel performance curves (cf. Eqs. 5a and 5b) needs to be showcased.
To achieve this objective, the OpenCPN model was utilised. This model can calculate sailboat routes with or without factoring in currents and incorporates shoreline knowledge, though it does not consider bathymetry. For our tests, we provided VISIR2 with wind and sea current fields identical to those used by OpenCPN (further details are provided in Sect. 5.1). Additionally, both models were equipped with the same sailboat polars. However, it is worth noting that OpenCPN does not handle leeway, whereas VISIR2 can manage it. The VISIR2 routes were computed on graphs of variable mesh resolution $\mathrm{1}/\mathrm{\Delta}x$ and connectivity ν, keeping fixed the “path resolution” parameter ΔP, which was introduced in Mannarini et al. (2019b, Eq. 6). This condition ensures that the maximum edge length remains approximately constant as the (Δx,ν) parameters are varied. Exemplary results are depicted in Fig. 9, with corresponding metrics provided in Table 7.
VISIR2 routes exhibit topological similarity to OpenCPN routes, yet for upwind sailing, they require a larger amount of tacking (see Fig. 9a). This discrepancy arises from the limited angular resolution of the graph (cf. Table 1). In the absence of currents, this implies a longer sailing time for VISIR2 with respect to OpenCPN routes, ranging between 1.0 % and 3.4 %. However, as the graph resolution is increased, the route duration decreases. Notably, this reduction plateaus at ν=7, indicating that such a resolution is optimal for the given path length. This type of comparison addresses the suggestion by Zis et al. (2020) to investigate the role of resolution in weather routing models. For a more thorough discussion of this particular aspect, please refer to Mannarini et al. (2019b).
Downwind routes all divert northwards because of stronger wind there (Fig. 9b). For these routes, the angular resolution does not pose a limiting factor, and VISIR2 routes exhibit shorter durations compared to OpenCPN routes. Considering the influence of currents as well, VISIR2 routes consistently prove to be the faster option, even for upwind sailing.
The disparities in duration between OpenCPN and VISIR2 routes could be attributed to various factors, including the interpolation method used for the wind field in both space and time and the approach employed to consider currents. Delving into these aspects would necessitate a dedicated investigation, which is beyond the scope of this paper.
Numerical tests have been integrated into the current VISIR2 release (Salinas et al., 2024a), covering the experiments listed in Tables 5–7 and beyond. These tests can be run using the Validazioni
module.
The computational performance of VISIR2 was evaluated using tests conducted on a single node of the “Juno” HPC facility at the EuroMediterranean Center on Climate Change (CMCC). This node was equipped with an Intel Xeon Platinum 8360Y processor, featuring 36 cores, each operating at a clock speed of 2.4 GHz, and boasting a pernode memory of 512 GB. Notably, parallelisation of the cores was not employed for these specific numerical experiments. Subsequently, our discussion here is narrowed down to assessing the performance of the module dedicated to computing optimal routes (“Tracce
”) in its motor vessel version.
In Fig. 10, we assess different variants of the shortestpath algorithm: leastdistance, leasttime, and leastCO_{2} procedures. We differentiate between the core of these procedures, which focuses solely on computing the optimal sequence of graph nodes (referred to hereafter as the “Dijkstra” component), and the broader procedure (“total”), which also includes the computation of both marine and vessel dynamical information along the legs of the optimal paths. The spatial interpolation option used in these tests (Sint =1 of Sect. 2.3.2) provides a conservative estimation of computational performance.
The numerical tests utilise the number of degrees of freedom (DOF) for the shortestpath problem as the independent variable. This value is computed as A⋅N_{τ}, where A denotes the number of edges and N_{τ} stands for the number of time steps of the fine grid (cf. Fig. 4). In the context of a seaonly edge graph, particularly in the case of a large graph (where border effects can safely be neglected), A can be represented as 4⋅N_{q1}(ν), where N_{q1} is defined by Eq. (20). Random edge weights were generated for graphs with ν=10, resulting in a number of DOFs ranging between 10^{5} and 10^{9}. Each data point in Fig. 10 represents the average of three identical runs, which helps reduce the impact of fluctuating HPC usage by other users. Additionally, the computational performance of VISIR2 is compared to that of VISIR1b, as documented in Mannarini and Carelli (2019b, Table 3, “With Tinterp”).
The primary finding is a confirmation of a powerlaw performance for all three optimisation objectives of VISIR2: distance, route duration, and total CO_{2} emissions. Remarkably, the curves appear nearly linear for the latter two algorithms (see Table 8). Such a scaling is even better than the linear–logarithmic worstcase estimate for Dijkstra's algorithms (Bertsekas, 1998, Sect. 2.3.1). Furthermore, this is not limited to just the Dijkstra components (as observed in VISIR1.b) but extends to the entire procedure, encompassing the reconstruction of alongroute variables. In addition to this enhanced scaling, VISIR2 demonstrates an improved absolute computational performance within the explored DOF range. The performance gain is approximately a factor of 10 when compared to VISIR1.b. This factor should even be larger when using the Sint =0 interpolation option.
Digging deeper into the details, we observe that the leastdistance procedure within VISIR2, while exclusively dealing with static edge weights, exhibits a less favourable scaling behaviour compared to both the leasttime and leastCO_{2} procedures. This is attributed to the postprocessing phase, wherein alongroute information has to be evaluated at the appropriate time step. Further development is needed to improve on this. Additionally, tiny nonlinearities are seen for the smaller number of DOFs. However, as proven by the metrics reported in Table 8, they do not affect the overall goodness of the powerlaw fits.
Lastly, it was found that peak memory allocation scales linearly across the entire explored range, averaging about 420 B per DOF. This is about 5 times larger than in VISIR1b and should be attributed to the NetworkX
structures used for graph representation. However, the large memory availability at the HPC facility prevented a possible degradation of performance for the largest numerical experiments due to memory swapping. A reduction in the unit memory allocation by a factor of 2 should be feasible using singleprecision floatingpoint format. Another strategy would involve using alternative graph libraries, such as igraph
.
A more comprehensive outcome of the VISIR2 code profiling, distinguishing also between the sailboat and the motor vessel version of the Tracce
module, is provided in Sect. S1 of the Supplement.
A prior version of VISIR2 has empowered GUTTAVISIR operational service, generating several million optimal routes within the Adriatic Sea over the span of a couple of years. In this section, we delve into outcomes stemming from deploying VISIR2 in different European seas. While the environmental fields are elaborated upon in Sect. 5.1, the results are given in Sect. 5.2, distinguishing by ferry and sailboat.
5.1 Environmental fields
The fields used for the case studies include both static and dynamic fields. The only static one was the bathymetry, extracted from the EMODnet product of 2022 (https://emodnet.ec.europa.eu/en/bathymetry, last access: 20 May 2024). Its spatial resolution was $\mathrm{1}/\mathrm{16}$ arcmin. The dynamic fields were the metocean conditions from both the European Centre for MediumRange Weather Forecasts (ECMWF) and CMEMS. Analysis fields from the ECMWF highresolution Atmospheric Model, 10 d forecast (Set I – HRES) with 0.1° resolution (https://www.ecmwf.int/en/forecasts/datasets/seti#Iia_fc, last access: 20 May 2024) were obtained. Both the u10m and v10m variables with 6hourly resolution were used. From CMEMS, analyses of the sea state, corresponding to the hourly MEDSEA_ANALYSISFORECAST_WAV_006_017 product, and of the sea surface circulation, corresponding to hourly MEDSEA_ANALYSISFORECAST_PHY_006_013, both with $\mathrm{1}/\mathrm{24}$° spatial resolution, were obtained. The windwave fields (vhm0_ww, vmdr_ww) and the Cartesian components (uo, vo) of the sea surface currents were used, respectively. Just for the comparison of VISIR2 to OpenCPN (see Sect. 3), 3hourly forecast fields from ECMWF (0.4°, https://www.ecmwf.int/en/forecasts/datasets/opendata, last access: 20 May 2024) and 3hourly RealTime Ocean Forecast System (RTOFS) forecasts ($\mathrm{1}/\mathrm{12}$°, https://polar.ncep.noaa.gov/global/about/, last access: 20 May 2024) for surface currents (p3049 and p3050 variables for U and V, respectively) were used.
The kinematics of VISIR2 presented in Sect. 2.1 do not limit the use to just surface ocean currents. This was just an initial approximation based on the literature discussed in Mannarini and Carelli (2019b). However, recent multisensor observations reported in Laxague et al. (2018) at a specific location in the Gulf of Mexico revealed a significant vertical shear, in both magnitude (by a factor of 2) and direction (by about 90°), within the first 8 m. Numerical ocean models typically resolve this layer; for instance the aforementioned Mediterranean product of CMEMS provides four levels within that depth. These vertically resolved data hold the potential to refine the computation of a ship's advection by the ocean flow. A plausible approach could involve the linear superposition of vessel velocity with a weighted average of the current, also considering the ship's hull geometry.
5.2 Results
To showcase some of the novel features of VISIR2, we present the outcomes of numerical experiments for both a ferry (as outlined in Sect. 2.5.1) and a sailboat (Sect. 2.5.2). All the results were generated using the interpolation options Sint =0 (as elaborated upon in Sect. 2.3.2) and Tint =2 (Sect. 2.3.1). These experiments considered the marine and atmospheric conditions prevailing in the Mediterranean Sea during the year 2022. Departures were scheduled daily at 03:00 UTC. The percentage savings (dQ) of a given quantity Q (such as the total CO_{2} emissions throughout the journey or the duration of sailing) are computed comparing the optimal (“opt”) to the leastdistance route (“gdt”):
5.2.1 Ferry
The chosen domain lies at the border between the Provençal Basin and the Ligurian Sea. Its sea state is significantly influenced by the mistral, a cold northwesterly wind that eventually affects much of the Mediterranean region during the winter months. The circulation within the domain is characterised by the southwestbound Liguro–Provençal current and the associated eddies (Schroeder and Chiggiato, 2022).
We conducted numerical experiments using VISIR2 with a graph resolution given by $(\mathit{\nu},\mathrm{1}/\mathrm{\Delta}x)=(\mathrm{4},\mathrm{12}/\mathit{\xb0})$, resulting in 2768 nodes and 114 836 edges within the selected domain. The time grid resolution was set at Δτ=30 min and N_{τ}=40. A single iteration (k=1) of Eq. (A1) was performed. The ferry engine load factor χ was varied to encompass values of 70 %, 80 %, 90 %, and 100 % of the installed engine power. For each day, both route orientations, with and without considering currents, were taken into account. This led to a total of 5840 numerical experiments. The computation time for each route was approximately 4 min, with the edge weight and shortestpath calculations consuming around 30 s.
In Fig. 11a an illustrative route is shown during a mistral event. As the ferry navigates against the wind, both its speed loss and its CO_{2} emission rate reach their maximum levels (cf. Fig. 6b, d). Consequently, both the leasttime and the leastCO_{2} algorithms calculate a detour into a calmer sea region where the combined benefits of improved sustained speed and reduced CO_{2} emissions compensate for the longer path's costs. The leastCO_{2} detour is wider than the leasttime one, as the additional duration is compensated for by greater CO_{2} savings. Moreover, these detours maximise the benefits from a southbound meander of the Liguro–Provençal current. Both optimal solutions intersect the water flow at points where it is narrowest, minimising the speed loss caused by the crosscurrent (cf. Eq. 16). Recessions of the isochrones become apparent starting at 6 h after departure. For this specific departure date and time, the overall reduction in CO_{2} emissions, in comparison to the shortestdistance route, exceeds 35 %.
Figure 11b illustrates that the magnitude of the related spatial diversion is merely intermediate, compared to the rest of 2022. Particularly during the winter months, the prevailing diversion is seen to occur towards the Ligurian Sea. Notably, VISIR2 even computed a diversion to the east of Corsica, which is documented in Sect. S3.1 of the Supplement. In the “Video supplement” accompanying this paper, all the 2022 routes between Porto Torres and Toulon are rendered, along with relevant environmental data fields.
To delve deeper into the statistical distribution of percentage CO_{2} savings defined as in Eq. (23), Fig. 12a provides a comparison with both the average significant wave height $\langle {H}_{\mathrm{s}}^{\left(\mathrm{gdt}\right)}\rangle $ and the absolute wave angle of attack $\langle \left{\mathit{\delta}}_{a}^{\left(\mathrm{gdt}\right)}\right\rangle $ along the shortestdistance route. Firstly, it should be noted that an increase in wave height can lead to either substantial or minimal CO_{2} emission savings. This outcome depends on whether the prevailing wave direction opposes or aligns with the vessel's heading. When focusing on routes with a CO_{2} saving of at least 2 %, it is seen that they mostly refer to either beam or head seas along the leastdistance route. This corresponds to elevated speed loss and subsequent higher emissions, as reported in Fig. 6b and d. This subset of routes shows a trend of larger savings in rougher sea states. Conversely, when encountering following seas with even higher H_{s}, savings remain below 1 %. This is due to both a smaller speed reduction and a lower CO_{2} emission rate. The counts of routes surpassing the 2 % saving threshold accounts for more than $\mathrm{1}/\mathrm{10}$ of the total routes, the ones above the 10 % threshold, represent about $\mathrm{1}/\mathrm{38}$th of the cases. This implies that, for the given ferry and the specified route, doubledigit percentage savings can be anticipated for about 10 calendar days per year.
The analysis of the CO_{2} savings distribution can be conducted by also considering the role of the engine load factor χ, as depicted in Fig. 12b. The distribution curves exhibit a biexponential shape, with the larger of the two decay lengths (d_{2}) inversely proportional to the magnitude of χ; cf. Table 10. This relationship is connected to the observation of reduced speed loss at higher χ as rougher sea conditions are experienced, which was already noted in the characteristics of this vessel in Sect. 2.5.1. The distribution's tail can extend to values ranging between 20 % and 50 %, depending on the specific value of χ.
Percentage CO_{2} savings, broken down by sailing direction and considering the presence or absence of currents, are detailed in Table 9. The average savings range from 0.6 % (1.0 % when considering sea currents) to 1.9 % (2.2 %). It is confirmed that the savings are more substantial on the route that navigates against the mistral wind (from Porto Torres to Toulon). However, the percentage savings are amplified when currents are taken into account, and this effect is particularly noticeable for the routes sailing in the downwind direction.
Further comments regarding the comparison of the CO_{2} savings presented here with the existing literature can be found in Sect. 6. While other works also hint at the presence of optimal route bundles, VISIR2 marks the first comprehensive exploration of how CO_{2} savings are distributed across various environmental conditions and engine loads.
5.2.2 Sailboat
The chosen area lies in the southern Aegean Sea, along a route connecting Greece (Monemvasia) and Türkiye (Marmaris). This area spans one of the most archipelagic zones in the Mediterranean Sea, holding historical significance as the birthplace of the term “archipelago”. The sea conditions in this area are influenced by the meltemi, prevailing northerly winds, particularly during the summer season. Such an “Etesian” weather pattern can extend its influence across a substantial portion of the Levantine basin (Lionello et al., 2008; Schroeder and Chiggiato, 2022). On the eastern side of the domain, the circulation is characterised by the westbound Asia Minor Current, while on its western flank, two prominent cyclonic structures separated by the west Cretan anticyclonic gyre are usually found (Theocharis et al., 1999).
We performed numerical experiments with VISIR2, with a graph resolution of $(\mathit{\nu},\mathrm{1}/\mathrm{\Delta}x)=(\mathrm{5},\mathrm{15}/\mathit{\xb0})$, leading to 2874 nodes and 156 162 edges in the selected domain. The resolution of the time grid was Δτ=30 min. Furthermore, N_{τ}=120 time steps of the environmental fields and k=2 iterations for Eq. (A1) were used. A First 36.7 sailboat was selected. For each day, both route orientations and all possible combinations of wind, current, and leeway were considered. This gave a total of 2920 numerical experiments. Each route required a total computing time of about 7 min, of which the edge weight and shortestpath computation amounted to 4 min, mainly spent on the edge weight computation. The excess time in comparison to the motor vessel's case study is attributed to both a higher value of N_{τ} and the additional time required for accounting for the exclusion of the nogo zone of the sailboat shown in Fig. 7.
In Fig. 13a, a sailboat route is depicted for a specific departure date, superimposed on the wind and sea current patterns. Both the leastdistance and leasttime routes appropriately steer clear of continental or insular landmasses, with the timeoptimal route opting for a more extensive detour. This adjustment is aimed at harnessing more favourable winds and circumventing unfavourable or cross currents, culminating in a remarkable 14.6 % reduction in route duration.
Moving to Fig. 13b, the collective set or “bundle” of eastbound routes is presented. Unlike the ferry routes showcased in Fig. 11b, it proves more challenging to discern a distinct seasonal pattern for the diversions of the sailboat routes. A metric for measuring diversions, such as the Fréchet distance utilised in Mannarini et al. (2019a), could facilitate the identification of patterns. The corresponding return routes are shown in Sect. S4.2 of the Supplement, confirming this trend. The bundles indicate that also accounting for currents leads to a more expansive set of optimal routes.
In only five cases (constituting $\mathrm{1.7}\times {\mathrm{10}}^{\mathrm{3}}$ of all sailboat routes), the leasttime route was discovered to be slower than the leastdistance route. These instances are scrutinised in Sect. S3.2 of the Supplement. The apparent inconsistency arises from the fact that these leasttime routes arrive at certain intermediate waypoints earlier but encounter less favourable sailing conditions compared to those encountered by the leastdistance routes arriving later. This discrepancy points to a nonFIFO situation (refer to Sect. 2.4.1). This scenario necessitates advancements in the leasttime algorithm to accommodate dynamic edge weights, potentially incorporating an initial waiting time, as discussed in Orda and Rom (1990).
A statistical evaluation of the time savings resulting from the optimisation process for sailboat routes is illustrated in Fig. 14a. Thereby, Eq. (23) is employed to assess both the path length and duration percentage savings. The −dT savings are generally proportional to path lengthening dL, particularly under nearly upwind conditions along the leastdistance route, i.e. where $\langle \left{\mathit{\delta}}_{i}^{\left(\mathrm{gdt}\right)}\right\rangle \sim {\mathit{\alpha}}_{\mathrm{0}}$ (cf. Eq. 7). This is understandable, as reduced sustained speeds and extended edge sailing times occur when wind originates from sectors close to the nogo zone, as depicted in Fig. 7a.
It is important to acknowledge that, under excessively weak or consistently sustained upwind conditions, a sailboat route might become unfeasible. A quantitative overview of such “failed” routes is provided in Table 11. It is evident that, thanks to the spatial diversions introduced by the route optimisation process, the likelihood of a leasttime route failing, compared to the leastdistance one, is reduced by a factor of approximately 100. In the “Video supplement” accompanying this paper, all sailboat routes between Monemvasia and Marmaris in 2022, along with the corresponding environmental fields, are included.
In Fig. 14b, the impact of currents and leeway is assessed. The influence of currents leads to a change in duration of up to approximately 5 % when compared to routes affected solely by the wind. Categorising the data based on sailing direction (as presented in Sect. S5 in the Supplement), currents primarily contribute to shorter route durations for westbound courses (benefiting from the Asia Minor Current). Conversely, they primarily result in extended durations for eastbound routes, particularly where, to the north of the island of Rhodes, there is no alternative but to sail against the current.
Turning to leeway, when not in combination with currents, it consistently extends the duration of routes, particularly, as indicated in the Supplement, when facing upwind conditions (more likely for westbound routes), as the speed loss is exacerbated due to a higher leeway speed (region with δ_{i}∼α_{0} in Fig. 7b).
As in our earlier comment in Sect. 2.1, the impact of leeway is mainly provided by its crosscourse component, which invariably decreases the vessel's SOG. Notably, the longitudinal component is smaller than the cross one by a tan δ factor; cf. Eq. (17). With δ estimated from Eq. (14) and Fig. 7b to fall within a range of a few degrees, the alongedge projection of leeway, ${w}_{\parallel}^{\left(L\right)}$, measures approximately $\mathrm{1}/\mathrm{10}$ of the transversal one, ${w}_{\perp}^{\left(L\right)}$.
When both effects, currents and leeway, are considered together, the distribution of duration changes in comparison to windonly routes resembles the distribution for the case with currents only. However, due to the impact of leeway, it is slightly skewed towards longer durations.
Finally in Table 11 time savings averaged throughout the year are presented. These savings are further categorised based on the direction of sailing and the specific combination of effects, including wind, currents, and leeway. The impact of sea currents generally increases the percentage duration savings from 2.4 % (the directional average of the windonly or wind and leeway cases) to 3.2 % (the average across all results affected by sea currents). We observe that routes featuring prevailing upwind conditions and favourable westbound currents, while also accounting for leeway, typically yield greater percentage duration savings compared to corresponding eastbound routes. This outcome can be attributed to the increase in duration of the leastdistance route, which results from the loss of sailboat manoeuvrability as the nogo zone is approached. It is also noted that the number of failed routes increases in the presence of leeway during upwind sailing. More statistical metrics are provided in Table S9 of the Supplement. Finally we observe that, in VISIR2, rigging is regarded as fixed (cf. Table 3), whereas in actual sailing practice, it may be optimised based on both the wind intensity and the point of sail.
In this section, we critically compare the CO_{2} savings achieved in Sect. 2.5.1 for the ferry's optimal routes with those reported in the ship weather routing literature (Sect. 6.1). Furthermore, we explore potential applications and users of VISIR2 (Sect. 6.2) and provide some insights into possible developments (Sect. 6.3).
6.1 Comparing CO_{2} savings with the literature
Only a handful of peerreviewed papers have reported results on emission savings through ship routing, yet a few of these findings are reviewed here for comparison with VISIR2's results.
For the ferry case study examined in this paper, the CO_{2} emissions can be halved compared to those along the leastdistance route in the best case (Fig. 12). This is in numerical agreement with the broad (0.1–48) % range reported in Bouman et al. (2017, Table 2). Notably, the upper limit reduction stands as an outlier, with the more probable values depicted in Bouman et al. (2017, Fig. 2) falling below 10 %. This aligns well with the results obtained from the statistical distribution of VISIR2 experiments, presented in both Fig. 12b and Table 9 of this paper.
Applying the VOIDS model, Mason et al. (2023a, Sect. 3.2) discovered that, for eastbound routes of a Panamax bulk carrier in the North Atlantic, voyage optimisation contributed to carbon savings ranging from 2.2 % to 13.2 %. This is a narrower range compared to the present findings of VISIR2 (cf. Fig. 12). However, both a different vessel type (a ferry) and a different domain of sailing (the western Mediterranean Sea) were considered in our numerical experiments.
Miola et al. (2011) presented data from the second IMO GHG study, where the estimated CO_{2} abatement potential for weather routing on the emissions of 2020 was reported to be as low as 0.24 %. In the same paper, a DNV study on projected emissions in 2030 was also cited, providing an estimate of 3.9 % for the CO_{2} abatement potential through weather routing. The former figure compares well to the average emission reduction computed via VISIR2 for the ferry downwind conditions and high engine load and the latter to results for upwind and low engine load (cf. Table 9).
Lindstad et al. (2013) estimated the reduction in CO_{2} emissions for a dry bulk Panamax vessel navigating in head seas during a typical stormy period in the North Atlantic. This reduction was determined when sailing on a 4500 nmi (nautical miles; 8300 km) route compared to the shorter (yet stormier) leastdistance route of 3600 nmi (6700 km). They found reductions ranging from 11 % to 48 %, depending on the speed of the vessel.
We note that, as for example in Mason et al. (2023a, Fig. 7), VISIR2 optimal routes also exhibit spatial variations contingent on the departure date, forming a bundle as illustrated in Fig. 11b and Sect. S4 of the Supplement. The shape of route boundaries was assessed for the United States’ Atlantic coast routes by means of AIS data in Breithaupt et al. (2017). However, while multimodal distributions depending on sailing direction were noted, the authors did not attribute the preferential lanes to the presence of ocean currents but speculated that it was due to bathymetric constraints or artificial aids to navigation.
VISIR possesses a capability to incorporate ocean currents into the voyage optimisation process. As shown in Mannarini and Carelli (2019b), this integration has proven to significantly reduce the duration of transatlantic routes. In the present paper, we reaffirm the positive impact of currents on ship route optimisation, also extending their benefits to the reduction in CO_{2} emissions (Table 9) and to the determination of more faithful duration savings for sailboat routes (Table 11).
In general, both average and extreme CO_{2} emission percentage savings found in the literature align well with the results obtained in the ferry case study presented in our paper. Nevertheless, engaging in a meaningful discussion of numerical differences, given the diverse range of vessel types, routes, environmental fields, and computational methods employed in the various published case studies, proves challenging.
VISIR2 contributes to the existing body of literature by providing an open computational platform that facilitates the simulation of optimal ship routes in the presence of waves, currents, and wind. These simulations are designed to be transparent, with customisable sea domain and vessel performance curves, allowing for thorough inspection, modification, and evaluation. This addresses the concern raised by Zis et al. (2020) regarding the necessity of benchmarking instances of optimal routes and the associated input data. By providing such benchmarks, VISIR2 supports and streamlines the work of future researchers in the field. Hence, we believe that the critical task of evaluating intermodel differences will best be addressed through dedicated intercomparison studies, as previously demonstrated with VISIR1 (Mannarini et al., 2019b).
6.2 Potential uses of VISIR2
Given its opensource nature, validated results, and numerical stability, VISIR2 can have great utility across various fields. The fact that both motor vessels and sailboats are treated equally will make VISIR2 suitable for use in weather routing of vessels with windassisted propulsion.
Moreover, VISIR2 can serve as a valuable tool for regulatory bodies seeking to make informed policies on shipping. As previously outlined, agencies like the European Maritime Safety Agency (EMSA), which oversee systems for monitoring, reporting, and verifying emissions, could utilise the model – provided that vessel performance curves and GHG emission rates are available – to calculate baseline emissions for various vessel and GHG types. With the inclusion of shipping in the EUETS, shipowners, ship managers, or bareboat charterers – whoever bears the fuel cost – are mandated to surrender allowances for their emissions. These stakeholders may find it beneficial to explore opensource solutions alongside existing commercial options.
VISIR2 also has the potential to help reduce uncertainty regarding the effectiveness of weather routing in reducing CO_{2} emissions (Bullock et al., 2020). For instance, evaluating distributions as in Fig. 12b, it becomes possible to characterise the joint potential of sea domains and vessel types for GHG emission savings. This concept also aligns with the idea of a “green corridor of shipping”, as envisioned by both the Clydebank Declaration (gov.uk, 2021) and the United States Department of State (DoS, 2022). In these initiatives, VISIR2, thanks to the generality of its optimisation algorithm (Algorithms 1, 2), could play a crucial role in minimising the consumption of costly zerocarbon fuel.
Furthermore, VISIR2 could be used to generate a dataset of optimal routes for the training of artificial intelligence systems for autonomous vessels (Li and Yang, 2023), surpassing the shortcomings of using AIS tracks, which include incomplete coverage (Filipiak et al., 2020). Finally, we note that, as opensource software, VISIR2 can even have educational purposes, providing training opportunities for ship officials and maritime surveillance authorities, as well as for beginner sailors.
6.3 Outlook
There are several possible avenues for future developments of VISIR2: the computer science, the algorithms, the ocean engineering, and the environmental fields.
First, as mentioned in Sect. 4, some computational performance improvements for the leastdistance procedure should be feasible. In applications where large domains, hyperresolution, or multipleinput environmental fields are required, it will be necessary to devise a solution that effectively reduces the computer's memory allocation. To further enhance modularity of VISIR2, future developments can focus on objectoriented programming principles. Containerisation of VISIR2 is currently underway as part of the development of a digital twin of the ocean (https://www.editomodellab.eu/, last access: 20 May 2024).
Further algorithmic work could address, for instance, construction of the graph directly in the projection space, facilitating the presence of more isotropic ship courses and perfectly collinear edges (cf. Sect. 2.2.2); development of an algorithm for givenduration, leastCO_{2} routes; incorporation of multiobjective optimisation techniques (Szlapczynska, 2015; Sidoti et al., 2017); generalisation of the leasttime algorithm to nonFIFO situations (cf. Sect. 2.4.1); and consideration of tacking time and motor assistance for sailboats or WAPSs.
Transitioning to upgrades in ocean engineering for VISIR2, the focus could shift towards targeting large oceangoing vessels, contingent upon the availability of related performance curves. Safety constraints on vessel intact stability (IMO, 2018b); considerations of slamming, green water, and lateral acceleration (Vettor and Guedes Soares, 2016); passenger comfort as highlighted by Carchen et al. (2021); or vessel performance in the presence of cross seas could be integrated. Secondgeneration intact stability criteria (Begovic et al., 2023) could be accounted for through a dynamic masking of the graph. VISIR2's readiness for windassisted ship propulsion hinges on the availability of an appropriate vessel performance curve that accounts for the added resistance from both wind and waves.
In terms of environmental data, it should be feasible to extend beyond the reliance on surface currents alone by incorporating the initial few depth layers of ocean models. Furthermore, VISIR2 currently operates under the assumption of having perfect knowledge of meteooceanographic conditions, which is provided through forecast fields for shorter voyages or analysis fields for longer ones. The latter corresponds to retracked routes as discussed in Mason et al. (2023b). However, for realtime applications during extended voyages, it is essential to incorporate adaptive routing strategies. This entails using the latest forecasts to execute rerouting as needed.
This paper presented the development of VISIR2: a modular, validated, and portable Pythoncoded model for ship weather routing. It provides a consistent framework for both motor vessels and sailboats by accounting for dynamic environmental fields such as waves, currents, and wind. The model can compute optimal ship routes even in complex and archipelagic domains. A cartographic projection, a feature that had been overlooked up to now, has been introduced in VISIR2. The model provides, for vessels with an angledependent performance curve, an improved level of accuracy in the velocity composition with sea currents. It is found that heading and course differ by an angle of attack, which is given by the solution of a transcendental equation (Eq. 13) involving an effective flow being the vector sum of currents and leeway. A computationally inexpensive iterative solution has been devised (Appendix A). Furthermore, a variant of Dijkstra's algorithm is introduced and used, which can minimise not just the CO_{2} emissions but any figure of merit depending on dynamic edge weights (cf. Algorithms 1, 2).
The validation of VISIR2 included comparisons to oracles and two intercomparison exercises. Differently from the few available ship weather routing packages or services, the VISIR2 software is accompanied by comprehensive documentation, making it suitable for community use.
The computational performance of the VISIR2 shortestpath module displayed a significant enhancement compared to its predecessor, VISIR1 (Sect. 4). A quasilinear scaling with problem complexity was demonstrated for up to 1×10^{9} DOFs. The robustness of VISIR2 was demonstrated across thousands of flawless route computations.
While the model is general with respect to vessel type, two case studies with VISIR2, based on realistic vessel seakeeping models, were documented in this paper.
From nearly 6000 routes of a 125 m long ferry, computed considering both waves and currents in the northwestern Mediterranean, average CO_{2} savings between 0.6 % and 2.2 %, depending on the engine load and prevailing sailing direction, were found. The distribution of the savings was biexponential, with the longer decay length becoming more pronounced at lower engine loads. This implied in particular that twodigit percentage CO_{2} savings, up to 49 %, were possible for about 10 d annually. This statistical distribution sheds new light on the underlying factors contributing to the variability observed in the role of weather routing, as reported in previous review studies (Bouman et al., 2017; Bullock et al., 2020). Furthermore, our findings bear significance for both the environmental impact of greenhouse gas emissions and the financial considerations within the EUETS.
From close to 3000 routes of an 11 m sailboat, within the southern Aegean Sea, accounting for both wind and currents, an average sailing time reduction of 2.4 % was observed. When considering currents as a factor, the duration of optimal routes could further be reduced by 3.2 %. Additionally, confirming prior work by Sidoti et al. (2023), disregarding the role of leeway would lead to an incorrect estimation of the route duration in upwind conditions. Several cases of nonFIFO behaviour were detected, and there is potential for addressing them in the future through the refinement of the current leasttime algorithm. All of these discoveries hold the potential to influence not just sailboat racing but also the utilisation of wind to aid in the propulsion of motor vessels.
In summary, this paper provides comprehensive documentation of the scientific hypotheses and decisions underpinning the development of an opensource ship routing model while also contributing to the quantification of achievable reductions in greenhouse gas emissions through voyage optimisation.
The angle δ between the ship's heading and course is obtained from the transcendental equation, Eq. (13). Its solution can be approximated by the following iteration:
where k is the number of iterations of the function
with γ being a constant resulting from the use of Eq. (2). The k=1 case corresponds to the solution provided in Mannarini and Carelli (2019b).
Both Eqs. (13) and (A1) were evaluated for a sailboat as in Sect. 2.5.2 using environmental conditions (wind, currents) for a domain in the central Adriatic Sea. A total of 11 hourly time steps and 18 474 edges from a graph with $(\mathit{\nu},\mathrm{1}/\mathrm{\Delta}x)=(\mathrm{4},\mathrm{12}/\mathit{\xb0})$ were considered, resulting in a total of about 2×10^{5} edge weight values. The iterative solution from Eq. (A1) was compared to the roots of Eq. (13) found via the scipy.optimize.root
solver, using as an initial guess the δ^{(1)} solution of Eq. (A1) (see the velocity_eval.py
function in the VISIR2 code). In what follows, the numerical solution from the solver is termed “exact”. The benefit of the approximated solution Eq. (A1) is that it can easily be parallelised on all graph arcs, while this is not possible for the exact solution, which processes one arc at a time.
The outcome for a sailboat is provided in Fig. A1a. It is seen that the iterative approximation departs from the exact solution for δ angles larger than about 5°. Such departures are mainly related to the effective cross flow ω_{⊥} (marker colour, determining the elongation from the origin). However, it is just a tiny fraction of the edges that present such departures, so the R^{2} correlation coefficient between the exact solution and its approximation is almost identical to 1 for any k>0, as shown in Fig. A1b.
The case k=0 corresponds to neglecting the loss of ships' momentum to balance the effective cross flow of Eq. (11b). Therefore, it wrongly underestimates the sailing times. For the First 36.7 sailboat under consideration here, the k=1 solution leads to a slope about 5 % off. Already for k=2 the correct slope is achieved within an error of 2 ‰. For the ferry, the k=1 iteration is sufficient to reach a 2 % accuracy; see Sect. S6.1 in the Supplement. This could be due to the ferry having a smoother angular dependence than the sailboat's one, as seen from Figs. 6b and 7a. Finally, in Fig. S22 of the Supplement, evidence of the validation of the exact solution in the absence of currents, Eq. (14), is also provided.
For identifying the vessel performance curves from a LUT via a neural network, a multilayer perceptron was used.
The models were built and trained via the scikitlearn
package (https://scikitlearn.org/stable/, last access: 20 May 2024). A 3fold crossvalidation was used to identify the best model for each vessel performance function. Different solvers, hidden layers' sizes, L2 regularisation terms, and activation functions were explored, covering a search space of about 10^{3} models. The optimal configuration made use of the rectified linear unit activation function, the Adam optimiser to minimise meansquared error, for at most 10^{3} passes through the training set (“epochs”) with a batch size of 200, a constant learning rate of 10^{−4}, and early stopping after the validation loss fails to decrease for 10 epochs.
The source code of VISIR2 is available from Salinas et al. (2024a) (https://doi.org/10.5281/zenodo.10960842). The distribution includes a user manual. Raw data in the form of input datasets and graphs used for the route computations are available from Salinas et al. (2024b) (https://doi.org/10.5281/zenodo.10674079). Intermediate data products in the form of routes for producing figures and tables in Sect. 5 are available from Salinas et al. (2024c) (https://doi.org/10.5281/zenodo.10674082).
Videos for this paper are available at https://doi.org/10.5446/s_1687 (ferry case study, Salinas, 2024a) and https://doi.org/10.5446/s_1688 (sailboat, Salinas, 2024b).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1743552024supplement.
GM: conceptualisation, funding acquisition, methodology, project administration, supervision, validation, writing (original draft), writing (review and editing); MLS: data curation, investigation, software, validation, visualisation; LC: data curation, investigation, software, validation, visualisation; NP: investigation, resources; JO: investigation, resources.
The contact author has declared that none of the authors has any competing interests.
The authors are not liable for casualties or losses that may occur in using routes computed via VISIR2 for navigation purposes.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Both ChatGPT3 and Gemini were utilised to review sections of this paper for Englishlanguage accuracy.
This research has been supported by Interreg (grant nos. 10043587 and 10253074) and Horizon Europe (grant nos. 101093293 and 101138583).
This paper was edited by David Ham and reviewed by two anonymous referees.
AlAboosi, F. Y., ElHalwagi, M. M., Moore, M., and Nielsen, R. B.: Renewable ammonia as an alternative fuel for the shipping industry, Current Opinion in Chemical Engineering, 31, 100670, https://doi.org/10.1016/j.coche.2021.100670, 2021. a
Begovic, E., Bertorello, C., Rinauro, B., and Rosano, G.: Simplified operational guidance for second generation intact stability criteria, Ocean Eng., 270, 113583, https://doi.org/10.1016/j.oceaneng.2022.113583, 2023. a
Bentley, J. L.: Multidimensional binary search trees used for associative searching, Communications of the ACM, 18, 509–517, 1975. a
Bertsekas, D.: Network Optimization: Continuous and Discrete Models, Athena Scientific, Belmont, Mass. 021789998, USA, 1998. a, b, c, d, e
Bouman, E. A., Lindstad, E., Rialland, A. I., and Strømman, A. H.: Stateoftheart technologies, measures, and potential for reducing GHG emissions from shipping – A review, Transport. Res. DTr. E., 52, 408–421, https://doi.org/10.1016/j.trd.2017.03.022, 2017. a, b, c, d
Breithaupt, S. A., Copping, A., Tagestad, J., and Whiting, J.: Maritime Route Delineation using AIS Data from the Atlantic Coast of the US, J. Navigation, 70, 379–394, https://doi.org/10.1017/S0373463316000606, 2017. a
Breivik, Ø. and Allen, A. A.: An operational search and rescue model for the Norwegian Sea and the North Sea, J. Marine Syst., 69, 99–113, 2008. a
Bullock, S., Mason, J., Broderick, J., and Larkin, A.: Shipping and the Paris climate agreement: a focus on committed emissions, BMC Energy, 2, 5, https://doi.org/10.1186/s42500020000152, 2020. a, b, c
Carchen, A., Gaggero, T., Besio, G., Mazzino, A., and Villa, D.: A method for the probabilistic assessment of the onboard comfort on a passenger vessel route, Ocean Eng., 225, 108702, https://doi.org/10.1016/j.oceaneng.2021.108, 2021. a
Claughton, A.: Developments in the IMS VPP Formulations, in: Fourteenth Chesapeake sailing yacht symposium, Annapolis, Maryland, 1–20, 1999. a
Claughton, A. R.: Developments in hydrodynamic force models for velocity prediction programs, in: Proceedings of the International Conference The Modern Yacht, The Royal Institution of Naval Architects, RINA, Paper: P20034 Proceedings, ISBN 0 903055 91 0, 2003. a
Dijkstra, E. W.: A note on two problems in connexion with graphs, Numerische Mathematik, 1.1, 269–271, https://doi.org/10.1145/3544585.3544600, 1959. a
DoS: Green Shipping Corridors Framework, Tech. rep., US Department of State, https://www.state.gov/greenshippingcorridorsframework/ (last access: 20 May 2024), 2022. a
Faber, J., van Seters, D., and Scholten, P.: Shipping GHG emissions 2030: Analysis of the maximum technical abatement potential, Tech. rep., CE Delft, 2023. a
Farkas, A., Parunov, J., and Katalinić, M.: Wave statistics for the middle Adriatic Sea, Pomorski zbornik, 52, 33–47, https://doi.org/10.18048/2016.52.02, 2016. a
Feeman, T. G.: Portraits of the Earth: A mathematician looks at maps, American Mathematical Soc., 18, 62–64, ISBN 0821832557, 2002. a
Filipiak, D., Węcel, K., Stróżyna, M., Michalak, M., and Abramowicz, W.: Extracting Maritime Traffic Networks from AIS Data Using Evolutionary Algorithm, Bus. Inf. Syst. Eng., 62, 435–450, https://doi.org/10.1007/s12599020006610, 2020. a
gov.uk: Clydebank Declaration, Tech. rep., UK Department for Transport, https://www.gov.uk/government/publications/cop26clydebankdeclarationforgreenshippingcorridors (last access: 20 May 2024), 2021. a
Guedes Soares, C.: Effect of heavy weather maneuvering on the waveinduced vertical bending moments in ship structures, J. Ship Res., 34, 60–68, 1990. a
IMO: MEPC.304(72) Initial IMO strategy on reduction of GHG emissions from ships, Tech. Rep. Annex 11, International Maritime Organization, London, UK, 2018a. a
IMO: SDC 5/J/7 Finalization of second generation intact stability criteria, Tech. rep., International Maritime Organization, London, UK, 2018b. a
IMO: MEPC.80/(WP.12) Report of the Working Group on Reduction of GHG Emissions from Ships Report of the Working Group on Reduction of GHG Emissions from Ships, Tech. rep., International Maritime Organization, London, UK, 2023. a
IPCC: Sixth Assessment Report, WG3, Ch.10, Tech. rep., IPCC, https://www.ipcc.ch/report/ar6/wg3/ (last access: 20 May 2024), 2022. a
IPCC: AR6 Synthesis Report: Climate Change 2023, Tech. rep., IPCC, https://www.ipcc.ch/report/ar6/syr/ (last access: 20 May 2024), 2023. a
Ladany, S. P. and Levi, O.: Search for optimal sailing policy, European J. Oper. Res., 260, 222–231, https://doi.org/10.1016/j.ejor.2016.12.013, 2017. a
Laxague, N. J., Özgökmen, T. M., Haus, B. K., Novelli, G., Shcherbina, A., Sutherland, P., Guigand, C. M., Lund, B., Mehta, S., Alday, M., and Molemaker, J.: Observations of nearsurface current shear help describe oceanic oil and plastic transport, Geophys. Res. Lett., 45, 245–249, 2018. a
Le Goff, C., Boussidi, B., Mironov, A., Guichoux, Y., Zhen, Y., Tandeo, P., Gueguen, S., and Chapron, B.: Monitoring the greater Agulhas Current with AIS data information, J. Geophys. Res.Oceans, 126, e2021JC017228, https://doi.org/10.1029/2021JC017228, 2021. a
Li, H. and Yang, Z.: Incorporation of AIS databased machine learning into unsupervised route planning for maritime autonomous surface ships, Transport. Res. ELog., 176, 103171, https://doi.org/10.1016/j.tre.2023.103171, 2023. a
Lindstad, H., Asbjørnslett, B. E., and Jullumstrø, E.: Assessment of profit, cost and emissions by varying speed as a function of sea conditions and freight market, Transport. Res. DTr. E., 19, 5–12, https://doi.org/10.1016/j.trd.2012.11.001, 2013. a
Lionello, P., Cogo, S., Galati, M., and Sanna, A.: The Mediterranean surface wave climate inferred from future scenario simulations, Global Planet. Change, 63, 152–162, https://doi.org/10.1016/j.gloplacha.2008.03.004, 2008. a
Lolla, S. V. T.: Path planning and adaptive sampling in the coastal ocean, Ph.D. thesis, Massachusetts Institute of Technology, http://hdl.handle.net/1721.1/103438 (last access: 20 May 2024), 2016. a
Maneewongvatana, S. and Mount, D. M.: It's okay to be skinny, if your friends are fat, in: Center for geometric computing 4th annual workshop on computational geometry, 2, 1–8, 1999. a
Mannarini, G. and Carelli, L.: [VISIR1.b ship routing model] source code (Matlab), Zenodo [code], https://doi.org/10.5281/zenodo.2563074, 2019a. a
Mannarini, G. and Carelli, L.: VISIR1.b: ocean surface gravity waves and currents for energyefficient navigation, Geosci. Model Dev., 12, 3449–3480, https://doi.org/10.5194/gmd1234492019, 2019b. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Mannarini, G., Lecci, R., and Coppini, G.: Introducing sailboats into ship routing system VISIR, in: 2015 6th International Conference on Information, Intelligence, Systems and Applications (IISA), IEEE, 1–6, https://doi.org/10.1109/IISA.2015.7387962, 2015. a
Mannarini, G., Pinardi, N., Coppini, G., Oddo, P., and Iafrati, A.: VISIRI: small vessels – leasttime nautical routes using wave forecasts, Geosci. Model Dev., 9, 1597–1625, https://doi.org/10.5194/gmd915972016, 2016a. a, b, c, d, e, f, g
Mannarini, G., Turrisi, G., D'Anca, A., Scalas, M., Pinardi, N., Coppini, G., Palermo, F., Carluccio, I., Scuro, M., Cretì, S., Lecci, R., Nassisi, P., and Tedesco, L.: VISIR: technological infrastructure of an operational service for safe and efficient navigation in the Mediterranean Sea, Nat. Hazards Earth Syst. Sci., 16, 1791–1806, https://doi.org/10.5194/nhess1617912016, 2016b. a
Mannarini, G., Carelli, L., Zissis, D., Spiliopoulos, G., and Chatzikokolakis, K.: Preliminary intercomparison of AIS data and optimal ship tracks, TransNav, 13, 53–61, https://doi.org/10.12716/1001.13.01.04, 2019a. a
Mannarini, G., Subramani, D., Lermusiaux, P., and Pinardi, N.: GraphSearch and Differential Equations for TimeOptimal Vessel Route Planning in Dynamic Ocean Waves, IEEE T. Intell. Transp., 21, 3581–3593, https://doi.org/10.1109/TITS.2019.2935614, 2019b. a, b, c, d, e, f
Mannarini, G., Carelli, L., Orović, J., Martinkus, C. P., and Coppini, G.: Towards LeastCO_{2} Ferry Routes in the Adriatic Sea, J. Marine Sci. Eng., 9, 115, https://doi.org/10.3390/jmse9020115, 2021. a
Mason, J., Larkin, A., Bullock, S., van der Kolk, N., and Broderick, J. F.: Quantifying voyage optimisation with wind propulsion for shortterm CO_{2} mitigation in shipping, Ocean Eng., 289, 116065, https://doi.org/10.1016/j.oceaneng.2023.116065, 2023a. a, b
Mason, J., Larkin, A., and GallegoSchmid, A.: Mitigating stochastic uncertainty from weather routing for ships with wind propulsion, Ocean Eng., 281, 114674, https://doi.org/10.1016/j.oceaneng.2023.114674, 2023b. a, b, c, d
Miola, A., Marra, M., and Ciuffo, B.: Designing a climate change policy for the international maritime transport sector: Marketbased measures and technological options for global and regional policy actions, Energy Policy, 39, 5490–5498, https://doi.org/10.1016/j.enpol.2011.05.013, 2011. a
Orda, A. and Rom, R.: Shortestpath and Minimumdelay Algorithms in Networks with Timedependent Edgelength, J. ACM, 37, 607–625, https://doi.org/10.1145/79147.214078, 1990. a, b, c
Salinas, M.: Ferry case study, TIB AVPortal [video], https://doi.org/10.5446/s_1687, 2024a. a
Salinas, M.: Sailboat case study, TIB AVPortal [video], https://doi.org/10.5446/s_1688, 2024b. a
Salinas, M. L., Carelli, L., and Mannarini, G.: [VISIR2 ship weather routing model] source code (Python), Zenodo [code], https://doi.org/10.5281/zenodo.10960842, 2024a. a, b, c, d, e
Salinas, M. L., Carelli, L., and Mannarini, G.: [VISIR2 ship weather routing model] raw data, Zenodo [data set], https://doi.org/10.5281/zenodo.10674079, 2024b. a
Salinas, M. L., Carelli, L., and Mannarini, G.: [VISIR2 ship weather routing model] intermediate products, Zenodo [data set], https://doi.org/10.5281/zenodo.10674082, 2024c. a
Schroeder, K. and Chiggiato, J.: Oceanography of the Mediterranean Sea: An Introductory Guide, Elsevier, ISBN 9780128236925, 2022. a, b
Sidoti, D., Avvari, G. V., Mishra, M., Zhang, L., Nadella, B. K., Peak, J. E., Hansen, J. A., and Pattipati, K. R.: A Multiobjective PathPlanning Algorithm With Time Windows for Asset Routing in a Dynamic WeatherImpacted Environment, IEEE T. Syst. Man Cyb., 47, 3256–3271, https://doi.org/10.1109/TSMC.2016.2573271, 2017. a
Sidoti, D., Pattipati, K. R., and BarShalom, Y.: Minimum Time Sailing Boat Path Algorithm, IEEE J. Ocean. Eng., 48, 307–322, https://doi.org/10.1109/JOE.2022.3227985, 2023. a, b
Smith, T. and Shaw, A.: An overview of the discussions from IMO MEPC 80 and Frequently Asked Questions, Tech. rep., UMAS, 2023. a
Svanberg, M., Ellis, J., Lundgren, J., and Landälv, I.: Renewable methanol as a fuel for the shipping industry, Renew. Sustain. Energ. Rev., 94, 1217–1228, https://doi.org/10.1016/j.rser.2018.06.058, 2018. a
Szlapczynska, J.: Multiobjective weather routing with customised criteria and constraints, J. Navigation, 68, 338–354, https://doi.org/10.1017/S0373463314000691, 2015. a
Tagliaferri, F., Philpott, A., Viola, I., and Flay, R.: On risk attitude and optimal yacht racing tactics, Ocean Eng., 90, 149–154, https://doi.org/10.1016/j.oceaneng.2014.07.020, 2014. a
Theocharis, A., Balopoulos, E., Kioroglou, S., Kontoyiannis, H., and Iona, A.: A synthesis of the circulation and hydrography of the South Aegean Sea and the Straits of the Cretan Arc (March 1994–January 1995), Prog. Oceanogr., 44, 469–509, https://doi.org/10.1016/S00796611(99)000415, 1999. a
van den Bremer, T. S. and Breivik, Ø.: Stokes drift, Philos. T. Roy. Soc. A, 376, 20170104, https://doi.org/10.1098/rsta.2017.0104, 2018. a
Vettor, R. and Guedes Soares, C.: Development of a ship weather routing system, Ocean Eng., 123, 1–14, https://doi.org/10.1016/j.oceaneng.2016.06.035, 2016. a, b
Wilson, G., Aruliah, D., Brown, C. T., Hong, N. P. C., Davis, M., Guy, R. T., Haddock, S. H., Huff, K. D., Mitchell, I. M., Plumbley, M. D., Waugh, B., White, E. P., and Wilson, P: Best practices for scientific computing, PLoS Biol., 12, e1001745, https://doi.org/10.1371/journal.pbio.1001745, 2014. a, b
Zis, T. P., Psaraftis, H. N., and Ding, L.: Ship weather routing: A taxonomy and survey, Ocean Eng., 213, 107697, https://doi.org/10.1016/j.oceaneng.2020.107697, 2020. a, b, c
 Abstract
 Introduction
 Technical advancements
 Validation
 Computational performance
 Case studies
 Discussion
 Conclusions
 Appendix A: Angle of attack
 Appendix B: Neural network features
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Article
(8420 KB)  Fulltext XML
 Companion paper 1
 Companion paper 2

Supplement
(12136 KB)  BibTeX
 EndNote
Ship weather routing has the potential to reduce CO_{2} emissions, but it currently lacks open and verifiable research. The Pythonrefactored VISIR2 model considers currents, waves, and wind to optimise routes. The model was validated, and its computational performance is quasilinear. For a ferry sailing in the Mediterranean Sea, VISIR2 yields the largest percentage emission savings for upwind navigation. Given the vessel performance curve, the model is generalisable across various vessel types.
Ship weather routing has the potential to reduce CO_{2} emissions, but it currently lacks open and...
 Abstract
 Introduction
 Technical advancements
 Validation
 Computational performance
 Case studies
 Discussion
 Conclusions
 Appendix A: Angle of attack
 Appendix B: Neural network features
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement