VISIR-2: ship weather routing in Python

Ship weather routing, which involves suggesting low-emission routes, holds potential for contributing to the decarbonisation of maritime transport. However, including because of a lack of readily deployable open-source and open-language computational models, its quantitative impact has been explored only to a limited extent.

As a response, the graph-search 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 K-dimensional tree. A least-CO₂ algorithm in the presence of dynamic graph edge weights has been implemented and validated, proving a quasi-linear 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, isochrone-bounded sectors have been introduced.

The resulting VISIR-2 model has been employed in numerical experiments within the Mediterranean Sea for the entirety of 2022, utilising meteo-oceanographic analysis fields. For a 125 m long ferry, the percentage saving of overall CO₂ expenditure follows a bi-exponential distribution. Routes with a carbon dioxide saving of at least 2 % with respect to the least-distance route were found for prevailing beam or head seas. Two-digit 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.

VISIR-2 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.

1 Introduction

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 mid-century, addressing all greenhouse gas (GHG) emissions, with a partial uptake of near-zero 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 (EU-ETS), starting from 2024 (https://climate.ec.europa.eu/eu-action/transport-emissions/reducing-emissions-shipping-sector_en, last access: 20 May 2024). For the first time ever, this will entail surrendering allowances for GHG emissions from vessels as well.

Zero-emission bunker fuels are projected to cost significantly more than present-day fossil fuels (Al-Aboosi 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. wind-assisted propulsion systems or WAPSs) and operational measures (e.g. speed reduction and ship weather routing). According to the “CE-Ship” 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 near-zero-GHG 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 (VISIR-1.a – Mannarini et al., 2016a, and VISIR-1.b – Mannarini and Carelli, 2019a). However, the use of a proprietary coding language (MATLAB) may hinder its further adoption. Also, the experience with VISIR-1 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 least-time routes were already estimated via VISIR-1.b, a dedicated algorithm for the direct computation of least-CO₂ routes was lacking.

To address all these requirements, we designed, coded, tested, and conducted extensive numerical experiments with the VISIR-2 model (Salinas et al., 2024a). VISIR-2 is a Python-coded software, inheriting from VISIR-1 the fact that it is based on a graph-search method. However, VISIR-2 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 VISIR-1. 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 least-CO₂ routes has been devised, building upon Dijkstra's original method. The visualisation of dynamic information along the route has been further improved using isochrone-bounded sectors. VISIR-2 features are thoroughly described in this paper, along with some case studies and suggestions for possible development lines in the future.

The VISIR-2 model enabled systematic assessment of CO₂ savings deriving from optimal use of meteo-oceanographic 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.

VISIR-2 is a numerical model well-suited 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/thetis-mrv.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 EU-ETS intensifies the importance of minimising GHG emissions. Internationally, even with the adoption of zero- or low-carbon fuels encouraged by the recent update of the IMO's strategy, it remains critical to save these fuels as much as possible. Therefore, VISIR-2 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 in-depth presentation of the innovations introduced by VISIR-2 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₂ emission percentage savings, some potential use cases, and an outlook on possible developments of VISIR-2 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 peer-reviewed 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 Real-time (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.avalon-routing.com/compute-route, 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 motor-assisted propulsion. Hourly departure dates within a couple of days are allowed, and ocean weather along the routes is provided in tabular form.

GUTTA-VISIR (https://gutta-visir.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 least-time and least-CO₂ 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 open-source platform, including a weather routing tool for sailboats. The product name originates from “Open-source 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 multi-objective (fuel consumption, duration, and safety) optimisation. An evolutionary algorithm was used, with initialisation from both random solutions and single-objective 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 dynamic-programming 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 dynamic-programming 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 dynamic-programming 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₂ savings has been provided. Furthermore, no related software has been disclosed in any case. The OpenCPN package lacks both a peer-reviewed publication and documentation of the methods implemented. A prevalence of dynamic-programming approaches is noted, especially for web applications, with the graph-search 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 open-source and well-documented 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.

2 Technical advancements

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 shortest-path 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 VISIR-2 code are presented in the software manual, provided along with its online repository (Salinas et al., 2024a).

2.1 Kinematics

For a graph-search 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 VISIR-2 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$\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.

Figure 1 Angular configuration with ($\mathit{\delta }=-\mathrm{27}$°, ${\mathit{\gamma }}_i=-\mathrm{109}$°), resulting in $\mathit{ϵ}=+\mathrm{1}$. The dark-grey area represents the ship hull, while the light-grey-shaded area denotes the no-go zone for α₀=25°. Clockwise-oriented arcs indicate positive angles, and filled circles at line ends denote the meteorological (“from”) convention.

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

$$\begin{array}{}\text{(1)}& \mathit{\delta }={\mathit{\psi }}_h-{\mathit{\psi }}_e,\end{array}$$

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

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}_f={\mathit{\psi }}_h-{\mathit{\psi }}_f=\mathit{\delta }-{\mathit{\gamma }}_f,\end{array}$$

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 VISIR-2 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

$$\begin{array}{}\text{(3)}& {\mathit{\psi }}_L={\mathit{\psi }}_h+\mathit{ϵ}\cdot \mathrm{90}\mathit{°},\end{array}$$

with

$$\begin{array}{}\text{(4)}& \mathit{ϵ}=\cos \left(\mathrm{180}\mathit{°}\lfloor {\mathit{\delta }}_i/\mathrm{180}\mathit{°}\rfloor \right),\end{array}$$

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

$$\begin{array}{}\text{(5a)}& {\displaystyle }F& {\displaystyle }=F\left(\right|{\mathit{\delta }}_i|,V_i;|{\mathit{\delta }}_a|,H_{\mathrm{s}},\mathit{\chi }),\text{(5b)}& {\displaystyle }L& {\displaystyle }=L\left(\right|{\mathit{\delta }}_i|,V_i),\end{array}$$

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

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_L=\mathrm{atan}(L/F).\end{array}$$

The above-introduced α_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 wind-magnitude-dependent no-go angle α₀(V_i), upwind navigation is not permitted, or

$$\begin{array}{}\text{(7)}& {\mathit{\delta }}_i\notin [-{\mathit{\alpha }}_{\mathrm{0}},{\mathit{\alpha }}_{\mathrm{0}}].\end{array}$$

Now, given a water flow expressed by the vector

$$\begin{array}{}\text{(8)}& \mathit{w}=C\widehat{w}=(u,v{)}^T\end{array}$$

and making reference to Fig. 1, the flow projections along ($\widehat{e}$) and across ($\widehat{o}$) the vehicle course are, respectively,

$$\begin{array}{}\text{(9a)}& {\displaystyle }& {\displaystyle }{w}_{\parallel }=C\cos ({\mathit{\psi }}_e-{\mathit{\psi }}_w)=u\sin \left({\mathit{\psi }}_e\right)+v\cos \left({\mathit{\psi }}_e\right),\text{(9b)}& {\displaystyle }& {\displaystyle }{w}_{\perp }=C\sin ({\mathit{\psi }}_e-{\mathit{\psi }}_w)=v\sin \left({\mathit{\psi }}_e\right)-u\cos \left({\mathit{\psi }}_e\right),\end{array}$$

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 cross-course projections of the leeway are given by

$$\begin{array}{}\text{(10a)}& {\displaystyle }{w}_{\parallel }^{\left(L\right)}& {\displaystyle }=L\cos ({\mathit{\psi }}_e-{\mathit{\psi }}_L)=-\mathit{ϵ}L\sin \mathit{\delta },\text{(10b)}& {\displaystyle }{w}_{\perp }^{\left(L\right)}& {\displaystyle }=L\sin ({\mathit{\psi }}_e-{\mathit{\psi }}_L)=-\mathit{ϵ}L\cos \mathit{\delta }.\end{array}$$

The simple relations on the right-hand side (r.h.s.) of Eqs. (10a) and (10b) follow from the similitude of the red- and green-shaded 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

$$\begin{array}{}\text{(11a)}& {\displaystyle }{\mathit{\omega }}_{\parallel }& {\displaystyle }=w_{\parallel }+w_{\parallel }^{\left(L\right)},\text{(11b)}& {\displaystyle }{\mathit{\omega }}_{\perp }& {\displaystyle }=w_{\perp }+w_{\perp }^{\left(L\right)}.\end{array}$$

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

$$\begin{array}{}\text{(12a)}& {\displaystyle }\text{SOG}=F\cos \left(\mathit{\delta }\right)+{\mathit{\omega }}_{\parallel },\text{(12b)}& {\displaystyle }\mathrm{0}=-F\sin \left(\mathit{\delta }\right)+{\mathit{\omega }}_{\perp }.\end{array}$$

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 angle-dependent 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:

$$\begin{array}{}\text{(13)}& \sin \mathit{\delta }={\displaystyle \frac{{\mathit{\omega }}_{\perp }(\mathit{\delta },{\mathit{\delta }}_i(\mathit{\delta }\left)\right)}{F\left(\right|{\mathit{\delta }}_i\left(\mathit{\delta }\right)\left)\right|,\left|{\mathit{\delta }}_a\right(\mathit{\delta }\left)\right|)}}\iff F\ne \mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(14)}& \mathit{\delta }=-\mathit{ϵ}{\mathit{\alpha }}_L.\end{array}$$

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 VISIR-2 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 VISIR-2.

We note that Eq. (13) holds if and only if

$$\begin{array}{}\text{(15)}& \left|{\mathit{\omega }}_{\perp }\right|\le F.\end{array}$$

Should this not be the case, the vessels' forward speed would not balance the effective drift.

As F is always non-negative, 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

$$\begin{array}{}\text{(16)}& \text{SOG}={\mathit{\omega }}_{\parallel }+\sqrt{F^{\mathrm{2}}-{\mathit{\omega }}_{\perp }^{\mathrm{2}}}.\end{array}$$

According to Eq. (16), the cross flow ω_⊥ always reduces SOG, as part of vehicle momentum must be spent balancing the drift. The along-edge 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 finite-difference quotient, the graph edge weight δt is computed as

$$\begin{array}{}\text{(17)}& \mathit{\delta }t={\displaystyle \frac{\mathit{\delta }x}{\text{SOG}}},\end{array}$$

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

$$\begin{array}{}\text{(18)}& \text{SOG}\ge \mathrm{0}\end{array}$$

should be checked in case the specific graph-search 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 counter-flow ω_∥ along a specific graph edge, whose weight must correspondingly be set to not a number.

The CO₂ emissions along a graph edge are given by

$$\begin{array}{}\text{(19)}& \mathit{\delta }{\mathrm{CO}}_{\mathrm{2}}=\mathrm{\Gamma }\mathit{\delta }t,\end{array}$$

where $\mathrm{\Gamma }=\mathrm{\Gamma }\left(\right|{\mathit{\delta }}_i|,V_i;|{\mathit{\delta }}_a|,H_{\mathrm{s}},\mathit{\chi })$ is the CO₂ emission rate of the specific vessel in the presence of the actual meteo-marine conditions. Both δt and δCO₂ are used in the least-CO₂ algorithm introduced in Sect. 2.4.

2.2 Graph generation

Graph preparation is crucial for any graph-search 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 graph-search method (such as Dijkstra's or A^*) and dynamic programming. Indeed, a dynamic-programming problem can be transformed into a shortest-path 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, VISIR-2's graph preparation incorporates several updates compared to its predecessor, VISIR-1b, as outlined in the following sub-sections.

Figure 2 Graph stencil for ν=4: (a) grid in spherical coordinates with Δx resolution along both latitude and longitude axes; (b) Mercator projection, with different resolutions along y or x, and graph edges (thick black and dashed grey lines) and angles (light blue) relative to due north. The y spacing is shown here as constant but, over a large latitudinal range, does vary. In VISIR-2, just the N_q1(ν) dark-grey nodes (cf. Eq. 20 and Table 1) are connected to the origin, while, in VISIR-1, all ν(ν+1) dark- or light-grey nodes were connected.

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, VISIR-2 introduces a projection feature, which was overlooked in both VISIR-1.a and VISIR-1.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 VISIR-2.

For a given sea domain, a graph is typically computed once and subsequently utilised for numerous different routes. However, in VISIR-1, the edge direction was recalculated every time a graph was utilised. In VISIR-2, 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 VISIR-2 graphs.

Finally it should be noted that in a directed graph, such as VISIR-2'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 Quasi-collinear 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 “quasi-collinear”. In VISIR-2, only the shortest one among these quasi-collinear 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

$$\begin{array}{}\text{(20)}& N_{q\mathrm{1}}\left(\mathit{\nu }\right)=\mathrm{2}\sum _{k=\mathrm{1}}^{\mathit{\nu }}\mathit{\phi }\left(k\right)\le \mathit{\nu }(\mathit{\nu }+\mathrm{1}),\end{array}$$

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 quasi-collinear ones. Using Eq. (20), at ν=4 already more than one-third of all edges are pruned and, at ν=10, nearly half of them are (cf. Table 1).

Table 1 Edge count and minimum angle with due north. The total number of edges in the first quadrant is ν(ν+1), and the count of non-collinear edges is N_q1 from Eq. (20). For the order of connectivity ν, the angle in the unprojected graph is given by $\mathrm{\Delta }\mathit{\theta }=\arcsin (\mathrm{1}/\mathit{\nu })$, regardless of latitude and grid resolution. Using a Mercator projection, at a latitude of L° and for Δx=D°, the angle is given by $\mathrm{\Delta }{\mathit{\theta }}_D^{\left(L\right)}$.

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 quasi-collinear 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.

Figure 3 Bathymetry field from EMODnet represented in shades of grey, with contour lines at depths at z=0 m and z=T, where T=7 m is the vessel draught. Additionally, the GSHHG shoreline at two different spatial resolutions is included.

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 VISIR-1b a procedure for restricting the check to inshore edges was introduced. In VISIR-2, as envisioned in Mannarini and Carelli (2019b, Appendix C), the process of searching for intersections is carried out using a K-dimensional 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 VISIR-2.

Various bathymetric databases can be used by VISIR-2. For European seas, the EMODnet dataset (https://portal.emodnet-bathymetry.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 “pseudo-shoreline” 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 pseudo-shoreline is the one used in VISIR-2 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 time-dependent environmental fields, whose information has to be remapped to the numerical grids of VISIR-2. This is done in a partly differently way than in VISIR-1, providing users with improved flexibility and control over numerical fields. These novel options are documented in what follows.

Figure 4 Temporal grid of VISIR-2. The upper horizontal axis (t) represents the coarse and uneven time resolution of the original environmental field. The lower horizontal axis corresponds to the fine and even time grid with resolution Δτ to which it is remapped.

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 VISIR-2, 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.

Figure 5 Spatial interpolation in VISIR-2. (a) The squares represent grid nodes of the environmental field φ(x), and the filled circles represent graph grid nodes. A graph edge is depicted as a magenta segment. (b) Transect of (a) along the edge direction $\widehat{e}$, with the interpolator of φ as a solid grey line. The (0, 1) subscripts refer to the value of the Sint parameter, while h and t refer to the edge head and tail, respectively.

VISIR-2 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 non-convex 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 VISIR-2.

The resulting edge-representative 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 so-called “sea-over-land” 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 Shortest-path algorithms

A major improvement made possible by the Python coding of VISIR-2 is the availability of built-in, advanced data structures such as dictionaries, queues, and heaps. They are key in the efficient implementations of graph-search algorithms (Bertsekas, 1998). In particular, as data structures are used, Dijkstra's algorithm worst-case performance can improve from quadratic, 𝒪(N²), 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

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathit{\delta }t\ge -\mathrm{1},\end{array}$$

an initial waiting time is not even unnecessary. This condition was assumed to hold in Mannarini et al. (2016a) for implementing a time-dependent Dijkstra's algorithm. That version of the shortest-path algorithm could not be used with an optimisation objective differing from voyage duration. As one aims to compute, for example, least-CO₂ routes, the algorithm requires further generalisation. This has been addressed in VISIR-2 via the pseudo-code 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 pseudo-function represents a NetworkX method to access the edge weight information.

We note the generality of the pseudo-code 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₂, 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 VISIR-2.

The shortest-distance and the least-time algorithms invoked for both motor vessels and sailboats are identical. Differences occur at the post-processing 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.

Algorithm 1 _DIJKSTRA_TDEP.

Require:

(G, source, target, wT, Ntau, Dtau), a NetworkX graph, source and target nodes, type of edge weight, maximum number of time steps, and time resolution, respectively

Ensure:

(costs, paths), Two dictionaries keyed by node ID: path costs from the source (e.g. cumulated CO₂) and corresponding optimal paths

1:

costs ←{}

2:

seen ←{ source:0 }

3:

paths ←{ source:[source] }

4:

# fringe is a min-priority queue of (cost, node) tuples

5:

fringe ← heap()

6:

fringe.push(0, source)

7:

while fringe ≠∅ do

8:

$(d,v)\leftarrow$ fringe.pop()

9:

if v∈ costs then

10:

# Already visited node

11:

skip

12:

end if

13:

costs[v] ←d

14:

if v=target and $\forall n\in G.\mathit{neigh}\left(\mathit{target}\right),n\in seen$ then

15:

exit

16:

end if

17:

$t\mathit{\_}idx\leftarrow \mathit{get}\mathit{\_}\mathit{time}\mathit{\_}\mathit{index}\left(\mathit{paths}\right[v],d,\mathit{wT},\mathit{Ntau},\mathit{Dtau})$

18:

# Iterate on v's forward-star

19:

for  (u, cost) in G.succ(v)  do

20:

# evaluate edge weight of wT type at time step t_idx

21:

$c\leftarrow \mathit{cost}.\mathit{at}\mathit{\_}\mathit{time}(t\mathit{\_}\mathit{idx},\mathit{wT})$

22:

$\mathit{vu}\mathit{\_}\mathit{cost}\leftarrow \mathit{costs}\left[v\right]+c$

23:

if  u∉seen or vu_cost<seen[u]  then

24:

seen[u]←vu_cost

25:

$\mathit{fringe}.\mathit{push}(\mathit{vu}\mathit{\_}\mathit{cost},u)$

26:

$\mathit{paths}\left[u\right]\leftarrow \mathit{paths}\left[v\right]+\left[u\right]$

27:

end if

28:

end for

29:

end while

Algorithm 2 GET_TIME_INDEX.

Require:

$(\mathit{paths},d,\mathit{wT},\mathit{Ntau},\mathit{Dtau})$, a dictionary of paths, node costs, type of edge weight, maximum number of time steps, and time resolution, respectively

Ensure:

t_idx, the time step at which the costs d are realised along the paths

1:

if wT= “time” then

2:

$t\mathit{\_}\mathit{idx}\leftarrow min(\mathit{Ntau},\lfloor d/\mathit{Dtau}\rfloor )$

3:

else

4:

# compute cTime cumulative time

5:

cTime←0

6:

t_idx←0

7:

for edge in paths do

8:

# evaluate edge delay at time step t_idx

9:

$\mathit{cTime}\leftarrow \mathit{cTime}+\mathit{edge}.\mathit{cost}.\mathit{at}\mathit{\_}\mathit{time}(t\mathit{\_}\mathit{idx},\mathrm{`}\mathrm{`}{\mathit{time}}^{\prime \prime })$

10:

$t\mathit{\_}\mathit{idx}\leftarrow \mathit{min}(\mathit{Ntau},\lfloor \mathit{time}/\mathit{Dtau}\rfloor )$

11:

end for

12:

end if

2.4.1 Non-FIFO

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 first-in–first-out (FIFO) strategy for traversing graph edges may not always be optimal. This is why such a scenario is referred to as “non-FIFO”. 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, non-FIFO edges were just 10⁻⁶ 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 no-go 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 least-time 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 VISIR-2 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 VISIR-1 the forward speed resulted, for motor vessels, from a semi-empirical 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 VISIR-2 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 VISIR-2 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 neural-network-based 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 VISIR-2 was a medium-size Ro-Pax vessel whose parameters are reported in Table 2.

Table 2 Principal parameters of the ferry.

1 kn = 1.852 km h⁻¹.

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

$$\begin{array}{}\text{(22)}& H_{\mathrm{s}}\left[\mathrm{m}\right]=\mathrm{0.0055}\cdot V_i\left[\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\right]+\mathrm{0.0127}\cdot {\left(V_i\left[\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\right]\right)}^{\mathrm{2}}.\end{array}$$

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 emission-related quantities. Leeway could not be considered by the simulator. The post-processed data feed the LUT to then be used for interpolating both STW and the CO₂ emission rate Γ as functions of significant wave height H_s, relative wind-wave 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₂ h⁻¹, shows a shallow dependence on H_s (Fig. 6c), while it is much more critically influenced by both χ and δ_a (Fig. 6c, d).

Figure 6 Ferry performance curve: in (a), STW is shown as a function of significant wave height H_s for head seas (δ_a=0°), with engine load χ indicated by the marker colour. In (b), STW is plotted as a function of δ_a at a constant χ=0.7, with H_s represented by the colour variation. The lower panels (c, d) display the CO₂ emission rate (Γ) with similar dependencies to in panels (a) and (b). Markers correspond to the LUT values, solid lines represent the spline interpolation, and dotted lines indicate the neural network's output.

2.5.2 Sailboat

Any sailboat described in terms of polars can in principle be used by VISIR-2. 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 four-degrees-of-freedom analysis, taking into account a wide range of semi-empirical hydrodynamic and aerodynamic models. It solves an equilibrium problem by a modified multi-dimensional Newton–Raphson iteration scheme. The analysis considered the added resistance due to waves by means of the so-called “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 wind-wave 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.

Table 3 Principal parameters of the sailboat (First 36.7).

Figure 7 Sailboat performance curve: forward speed F in (a) and leeway speed L in (b) are both plotted against the true wind angle δ_i. Panel (c) shows the leeway angle α_L obtained from Eq. (6). Marker and line colours represent wind magnitude V_i. Data start at δ_i=α₀(V_i). Markers refer to the LUT, solid lines to spline interpolation, and dotted lines to the neural network's output. The colour bar also reports the LUT's minimum and maximum values printed in blue and red, respectively.

The outcome corresponds to Eqs. (5a) and (5b) and is provided in Fig. 7. The no-go angle α₀ 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 no-go zone (δ_i=α₀) and beam reach (δ_i=90°); see Fig. 7b. As the point of sail transitions from the no-go 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 VISIR-2 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 stripe-wise rendering introduced in VISIR-1.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 VISIR-2.

In addition to isochrones, lines of equal distance from the route's origin (or “isometres”) and lines of equal quantities of CO₂ 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₂ 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 VISIR-2. While in VISIR-1 modularity was limited to the graph preparation, which was detached from the main pipeline (Mannarini et al., 2016a, Fig. 8), the VISIR-2 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.

Figure 8 VISIR-2 workflow. Modules enclosed within thicker frames are intended for direct execution by the end user, while the other modules can be modified for advanced usage. The data flow occurs along the wavy arrow, with routine calls along the dashed line.

Table 4 VISIR-2 modules with their original names, purpose, and references within this paper. Module nos. 1–5 represent the core package and require the visir-venv virtual environment. Module no. 6 runs with visir_vis-venv.

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 GUTTA-VISIR) instead, the computational workflow can be streamlined by avoiding the saving of the intermediate results. VISIR-2 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 VISIR-2 released here, it facilitates the ports' selection from the World Port Index (https://msi.nga.mil/Publications/WPI, last access: 20 May 2024) database.

VISIR-2 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 high-performance computing (HPC) facilities.

3 Validation

Validating a complex model like VISIR-2 is imperative. The code was developed with specific runs of VISIR-1 as a benchmark. The validation of VISIR-1 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 VISIR-1 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 VISIR-2, at different graph resolutions, as shown in Table 5. The errors are consistently found to be below 1 %.

Table 5 VISIR-2 route durations compared to analytic oracles (Cycloid and Techy), as referenced in the main text. L₀ and T₀ represent the length scales and timescales, respectively. The oracle durations are denoted by T^(e), while those from VISIR-2 are denoted by T^*. The percentage mismatch is calculated as $d{T}^{*}=(T^{*}/T^{\left(e\right)})-\mathrm{1}$.

1 nmi = 1852 m.

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 VISIR-2, the benchmarks were taken as LSE simulations at the nearest grid resolution. Notably, VISIR-2 consistently produces shorter-duration routes compared to VISIR-1.b (see Table 6).

Table 6 VISIR-2 vs. LSE durations. The relative error $\mathrm{d}{T}_{\mathrm{res}}^{*}$ is defined as the discrepancy between T^* and LSE at two different grid resolutions $\mathrm{1}/\mathrm{\Delta }x$. Both VISIR-2 and VISIR-1 outcomes are provided.

For both analytical and numerical benchmarks, distinct from the scenario discussed in Sect. 2.2.3, quasi-collinear 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 VISIR-2 to accommodate angle-dependent 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 VISIR-2 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 VISIR-2 can manage it. The VISIR-2 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.

Table 7 VISIR-2 vs. OpenCPN comparison. Durations T^* and relative mismatch dT^* for the cases shown in Fig. 9 are provided. k=2 and Δτ=15 min used throughout the numerical experiments.

Figure 9 VISIR-2 routes with wind and currents vs. OpenCPN: graphs of variable resolution, indexed by ν as shown in the legend, with a constant ΔP∼0.3°. Field intensity is in grey tones, and the direction is shown as black streamlines. Shell representation with isochrones is in dashed gold lines, and labels are in hours. The OpenCPN solution is plotted as a navy line. Panels (a) and (b) refer to the west- and eastbound voyage, respectively.

VISIR-2 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 VISIR-2 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 VISIR-2 routes exhibit shorter durations compared to OpenCPN routes. Considering the influence of currents as well, VISIR-2 routes consistently prove to be the faster option, even for upwind sailing.

The disparities in duration between OpenCPN and VISIR-2 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 VISIR-2 release (Salinas et al., 2024a), covering the experiments listed in Tables 5–7 and beyond. These tests can be run using the Validazioni module.

4 Computational performance

The computational performance of VISIR-2 was evaluated using tests conducted on a single node of the “Juno” HPC facility at the Euro-Mediterranean 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 per-node 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 shortest-path algorithm: least-distance, least-time, and least-CO₂ 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 shortest-path 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 sea-only 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⁵ and 10⁹. 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 VISIR-2 is compared to that of VISIR-1b, as documented in Mannarini and Carelli (2019b, Table 3, “With T-interp”).

Figure 10 Profiling of computing time for the Tracce module (motor vessel case). The independent variable is the number of DOFs (#DOF) in the graph. Markers refer to experimental data points and lines to least-square fits. Void markers and dashed lines refer to just Dijkstra's component, while full markers and solid lines refer to the whole of the each routine. The colours refer to the three alternative optimisation objectives, while black is used for VISIR-1.b results.

The primary finding is a confirmation of a power-law performance for all three optimisation objectives of VISIR-2: distance, route duration, and total CO₂ 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 worst-case estimate for Dijkstra's algorithms (Bertsekas, 1998, Sect. 2.3.1). Furthermore, this is not limited to just the Dijkstra components (as observed in VISIR-1.b) but extends to the entire procedure, encompassing the reconstruction of along-route variables. In addition to this enhanced scaling, VISIR-2 demonstrates an improved absolute computational performance within the explored DOF range. The performance gain is approximately a factor of 10 when compared to VISIR-1.b. This factor should even be larger when using the Sint =0 interpolation option.

Digging deeper into the details, we observe that the least-distance procedure within VISIR-2, while exclusively dealing with static edge weights, exhibits a less favourable scaling behaviour compared to both the least-time and least-CO₂ procedures. This is attributed to the post-processing phase, wherein along-route 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 power-law 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 VISIR-1b 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 single-precision floating-point format. Another strategy would involve using alternative graph libraries, such as igraph.

Table 8 Fit coefficients of the $T_c=a\cdot {\text{DOF}}^b+c$ regressions for various components of Tracce, motor vessel version. “D” stands for Dijkstra's algorithm only, while “tot” includes the post-processing for reconstructing the voyage. p(K) is the p value for the K coefficient. All data refer to VISIR-2 except the *_V1b data, which refer to VISIR-1.b.

A more comprehensive outcome of the VISIR-2 code profiling, distinguishing also between the sailboat and the motor vessel version of the Tracce module, is provided in Sect. S1 of the Supplement.

5 Case studies

A prior version of VISIR-2 has empowered GUTTA-VISIR 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 VISIR-2 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 Medium-Range Weather Forecasts (ECMWF) and CMEMS. Analysis fields from the ECMWF high-resolution Atmospheric Model, 10 d forecast (Set I – HRES) with 0.1° resolution (https://www.ecmwf.int/en/forecasts/datasets/set-i#I-i-a_fc, last access: 20 May 2024) were obtained. Both the u10m and v10m variables with 6-hourly 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 wind-wave fields (vhm0_ww, vmdr_ww) and the Cartesian components (uo, vo) of the sea surface currents were used, respectively. Just for the comparison of VISIR-2 to OpenCPN (see Sect. 3), 3-hourly forecast fields from ECMWF (0.4°, https://www.ecmwf.int/en/forecasts/datasets/open-data, last access: 20 May 2024) and 3-hourly Real-Time 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 VISIR-2 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 multi-sensor 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 VISIR-2, 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₂ emissions throughout the journey or the duration of sailing) are computed comparing the optimal (“opt”) to the least-distance route (“gdt”):

$$\begin{array}{}\text{(23)}& \mathrm{d}Q={\displaystyle \frac{Q^{\left(\mathrm{opt}\right)}-Q^{\left(\mathrm{gdt}\right)}}{Q^{\left(\mathrm{gdt}\right)}}}.\end{array}$$

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 southwest-bound Liguro–Provençal current and the associated eddies (Schroeder and Chiggiato, 2022).

We conducted numerical experiments using VISIR-2 with a graph resolution given by $(\mathit{\nu },\mathrm{1}/\mathrm{\Delta }x)=(\mathrm{4},\mathrm{12}/\mathit{°})$, 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 shortest-path 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₂ emission rate reach their maximum levels (cf. Fig. 6b, d). Consequently, both the least-time and the least-CO₂ algorithms calculate a detour into a calmer sea region where the combined benefits of improved sustained speed and reduced CO₂ emissions compensate for the longer path's costs. The least-CO₂ detour is wider than the least-time one, as the additional duration is compensated for by greater CO₂ 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₂ emissions, in comparison to the shortest-distance 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, VISIR-2 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.

Figure 11 The ferry's optimal routes between Porto Torres (ITPTO) and Toulon (FRTLN) with both waves and currents: (a) for the specified departure date and time, the shortest-distance route is shown in blue, the least-time route in red, and the least-CO₂ route in green. The H_s field is displayed in shades of grey with black arrows, while the currents are depicted in purple tones with white streamlines. The shortest-path algorithm did not utilise environmental field values within the etched area. Additionally, isochrones of the CO₂-optimal route are shown at 3-hourly intervals. The engine load was χ=0.7. (b) A bundle of all northbound CO₂-optimal routes (for $\mathit{\chi }=\mathit{\{}\mathrm{0.7},\mathrm{0.8},\mathrm{0.9},\mathrm{1.0}\mathit{\}}$) is presented, with the line colour indicating the departure month.

To delve deeper into the statistical distribution of percentage CO₂ 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 shortest-distance route. Firstly, it should be noted that an increase in wave height can lead to either substantial or minimal CO₂ 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₂ saving of at least 2 %, it is seen that they mostly refer to either beam or head seas along the least-distance 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₂ 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, double-digit percentage savings can be anticipated for about 10 calendar days per year.

The analysis of the CO₂ 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 bi-exponential shape, with the larger of the two decay lengths (d₂) 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 χ.

Table 9 Average relative savings of the CO₂-optimal vs. the least-distance route (in %), for various engine loads (χ), considering just waves (wa) or also currents (wa-cu), for ferry routes between Toulon (FRTLN) and Porto Torres (ITPTO) as in Fig. 11. The χ-averaged values are also provided in the “Avg” columns.

Percentage CO₂ 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.

Figure 12 Metrics relative to ferry routes pooled on sailing directions (FRTLN ↔ ITPTO) and χ, using both waves and currents. (a) Percentage savings, with the markers' grey shade representing the mean angle of attack along the least-distance route. The total number of routes, including those with relative CO₂ savings above 2 % (solid line) and 10 % (dashed), are also provided; (b) Distributions of the CO₂ savings for each χ value, with fitted bi-exponential functions as in Table 10. Each set of four columns pertains to a bin centred on the nearest tick mark and spanning a width of 5 %.

Table 10 Fit coefficients of $y=a\cdot \left[\exp \right(-x/d_{\mathrm{1}})+b\cdot \exp (-x/d_{\mathrm{2}}\left)\right]$ on the data of Fig. 12b.

Further comments regarding the comparison of the CO₂ 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, VISIR-2 marks the first comprehensive exploration of how CO₂ 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 VISIR-2, with a graph resolution of $(\mathit{\nu },\mathrm{1}/\mathrm{\Delta }x)=(\mathrm{5},\mathrm{15}/\mathit{°})$, 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 shortest-path 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 no-go 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 least-distance and least-time routes appropriately steer clear of continental or insular landmasses, with the time-optimal 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 least-time route was discovered to be slower than the least-distance route. These instances are scrutinised in Sect. S3.2 of the Supplement. The apparent inconsistency arises from the fact that these least-time routes arrive at certain intermediate waypoints earlier but encounter less favourable sailing conditions compared to those encountered by the least-distance routes arriving later. This discrepancy points to a non-FIFO situation (refer to Sect. 2.4.1). This scenario necessitates advancements in the least-time algorithm to accommodate dynamic edge weights, potentially incorporating an initial waiting time, as discussed in Orda and Rom (1990).

Figure 13 The sailboat's optimal routes between GRMON and TRMRM, considering both wind and currents: (a) for the specified departure date and time, the least-time route is depicted in red, and the shortest-distance one is in blue. The wind field is represented in shades of grey with black arrows, while the currents are shown in purple tones with white streamlines. Additionally, isochrones of the time-optimal route are displayed at 6-hourly intervals. (b) A bundle of all eastbound time-optimal routes is presented, with the line colour indicating the departure month.

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 least-distance 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 no-go 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 least-time route failing, compared to the least-distance 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∼α₀ in Fig. 7b).

As in our earlier comment in Sect. 2.1, the impact of leeway is mainly provided by its cross-course 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 along-edge 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 wind-only 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 wind-only 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 least-distance route, which results from the loss of sailboat manoeuvrability as the no-go 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 VISIR-2, 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.

Figure 14 Metrics of the sailboat's optimal routes pooled on sailing directions (GRMON ↔ TRMRM). (a) Duration percentage savings −dT vs. relative lengthening dL considering just wind: the markers' grey shade represents the average angle of attack of wind $|\langle {\mathit{\delta }}_i^{\left(\mathrm{gdt}\right)}\rangle |$ along the least-distance route. The number of data (#data) is given by ${\sum }_d\mathrm{365}-N_f^{\left(g\right)}$, where d represents the two sailing directions and the $N_f^{\left(g\right)}$ values are from the first row in Table 11. (b) Histograms of relative route duration T_f, with environmental forcing combination f defined by the column colour, with respect to the duration T_wi of the wind-only optimal routes.

Table 11 Average time savings of the sailboat routes (in %), considering just wind (wi) or also various combinations of currents (cu) and leeway (le) for the sailboat routes between Monemvasia (GRMON) and Marmaris (TRMRM) as in Fig. 13. The number of failed routes for the least-distance $N_f^{\left(g\right)}$ or the least-time routes $N_f^{\left(o\right)}$ is also provided.

6 Discussion

In this section, we critically compare the CO₂ 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 VISIR-2 (Sect. 6.2) and provide some insights into possible developments (Sect. 6.3).

6.1 Comparing CO₂ savings with the literature

Only a handful of peer-reviewed papers have reported results on emission savings through ship routing, yet a few of these findings are reviewed here for comparison with VISIR-2's results.

For the ferry case study examined in this paper, the CO₂ emissions can be halved compared to those along the least-distance 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 VISIR-2 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 VISIR-2 (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₂ 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₂ abatement potential through weather routing. The former figure compares well to the average emission reduction computed via VISIR-2 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₂ 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) least-distance 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), VISIR-2 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₂ emissions (Table 9) and to the determination of more faithful duration savings for sailboat routes (Table 11).

In general, both average and extreme CO₂ 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.

VISIR-2 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, VISIR-2 supports and streamlines the work of future researchers in the field. Hence, we believe that the critical task of evaluating inter-model differences will best be addressed through dedicated inter-comparison studies, as previously demonstrated with VISIR-1 (Mannarini et al., 2019b).

6.2 Potential uses of VISIR-2

Given its open-source nature, validated results, and numerical stability, VISIR-2 can have great utility across various fields. The fact that both motor vessels and sailboats are treated equally will make VISIR-2 suitable for use in weather routing of vessels with wind-assisted propulsion.

Moreover, VISIR-2 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 EU-ETS, 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 open-source solutions alongside existing commercial options.

VISIR-2 also has the potential to help reduce uncertainty regarding the effectiveness of weather routing in reducing CO₂ 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, VISIR-2, thanks to the generality of its optimisation algorithm (Algorithms 1, 2), could play a crucial role in minimising the consumption of costly zero-carbon fuel.

Furthermore, VISIR-2 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 open-source software, VISIR-2 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 VISIR-2: the computer science, the algorithms, the ocean engineering, and the environmental fields.

First, as mentioned in Sect. 4, some computational performance improvements for the least-distance procedure should be feasible. In applications where large domains, hyper-resolution, or multiple-input 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 VISIR-2, future developments can focus on object-oriented programming principles. Containerisation of VISIR-2 is currently underway as part of the development of a digital twin of the ocean (https://www.edito-modellab.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 given-duration, least-CO₂ routes; incorporation of multi-objective optimisation techniques (Szlapczynska, 2015; Sidoti et al., 2017); generalisation of the least-time algorithm to non-FIFO 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 VISIR-2, the focus could shift towards targeting large ocean-going 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. Second-generation intact stability criteria (Begovic et al., 2023) could be accounted for through a dynamic masking of the graph. VISIR-2's readiness for wind-assisted 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, VISIR-2 currently operates under the assumption of having perfect knowledge of meteo-oceanographic 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 real-time applications during extended voyages, it is essential to incorporate adaptive routing strategies. This entails using the latest forecasts to execute re-routing as needed.

7 Conclusions

This paper presented the development of VISIR-2: a modular, validated, and portable Python-coded 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 VISIR-2. The model provides, for vessels with an angle-dependent 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₂ emissions but any figure of merit depending on dynamic edge weights (cf. Algorithms 1, 2).

The validation of VISIR-2 included comparisons to oracles and two inter-comparison exercises. Differently from the few available ship weather routing packages or services, the VISIR-2 software is accompanied by comprehensive documentation, making it suitable for community use.

The computational performance of the VISIR-2 shortest-path module displayed a significant enhancement compared to its predecessor, VISIR-1 (Sect. 4). A quasi-linear scaling with problem complexity was demonstrated for up to 1×10⁹ DOFs. The robustness of VISIR-2 was demonstrated across thousands of flawless route computations.

While the model is general with respect to vessel type, two case studies with VISIR-2, 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₂ savings between 0.6 % and 2.2 %, depending on the engine load and prevailing sailing direction, were found. The distribution of the savings was bi-exponential, with the longer decay length becoming more pronounced at lower engine loads. This implied in particular that two-digit percentage CO₂ 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 EU-ETS.

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 non-FIFO behaviour were detected, and there is potential for addressing them in the future through the refinement of the current least-time 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 open-source ship routing model while also contributing to the quantification of achievable reductions in greenhouse gas emissions through voyage optimisation.

Appendix A: Angle of attack

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:

$$\begin{array}{}\text{(A1)}& \begin{array}{lll}{\mathit{\delta }}^{\left(\mathrm{0}\right)}& =& \mathrm{0},\\ {\mathit{\delta }}^{\left(k\right)}& =& h\left({\mathit{\delta }}^{(k-\mathrm{1})}\right)\mathrm{for}k=\mathrm{1},\mathrm{2},\mathrm{\dots },\end{array}\end{array}$$

where k is the number of iterations of the function

$$\begin{array}{}\text{(A2)}& h\left(x\right)={\displaystyle }\arcsin \left({\displaystyle \frac{{\mathit{\omega }}_{\perp }(\mathit{\delta }=x,{\mathit{\delta }}_i=x-\mathit{\gamma })}{F\left(\right|x-\mathit{\gamma }\left|\right)}}\right),\end{array}$$

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).

Figure A1 Approximate vs. exact solution of Eq. (13) for a First 36.7 sailboat. (a) Iterative solution of Eq. (A1) with k=1 vs. the exact solution, using ω_⊥ as the marker colour; (b) unexplained variance (R is the Pearson's correlation coefficient) of the linear regression and fitted slope coefficient for various k values.

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{°})$ were considered, resulting in a total of about 2×10⁵ 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 δ⁽¹⁾ solution of Eq. (A1) (see the velocity_eval.py function in the VISIR-2 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² 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.

Appendix B: Neural network features

For identifying the vessel performance curves from a LUT via a neural network, a multi-layer perceptron was used.

The models were built and trained via the scikit-learn package (https://scikit-learn.org/stable/, last access: 20 May 2024). A 3-fold cross-validation 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³ models. The optimal configuration made use of the rectified linear unit activation function, the Adam optimiser to minimise mean-squared error, for at most 10³ passes through the training set (“epochs”) with a batch size of 200, a constant learning rate of 10⁻⁴, and early stopping after the validation loss fails to decrease for 10 epochs.