VISIR-I . b : waves and ocean currents for energy efficient navigation

VISIR-I.b, the latest development of the ship routing model published in Mannarini et al. (2016a), is here presented. The new model version targets large ocean-going vessels by accounting for both waves and ocean currents. In order to effectively use currents in a graph-search method, new equations are derived and validated versus analytical benchmarks. A case study is computed in the Atlantic Ocean, on a route from the Chesapeake Bay to the Mediterranean Sea and vice versa. Ocean analysis fields from data-assimilative models (for both ocean state and hydrodynamics) are employed. The impact of 5 waves and ocean currents on transatlantic crossings is assessed through mapping of the spatial variability of the routes, analysis of their kinematics, distribution of the optimal voyage duration vs. its length, and impact on the Energy Efficiency Operational Indicator of the International Maritime Organization. It is distinguished between sailing with or against the main ocean current. The seasonal dependence of the savings is evaluated, indicating, for the featured case study, larger savings during the summer crossings and larger intra-monthly variability in winter. The monthly-mean savings sum up to values between 3 and 12%, 10 while the contribution of ocean currents is between 1 and 4%. Also, several other ocean routes are considered, providing a pan-Atlantic scenario assessment of the potential gains in energy efficiency from optimal tracks and linking them to regional meteo-oceanographic features.


Introduction
The strongest water flows are generally observed in ocean Western boundary currents, in tidal currents, in the circulation of straits and fjords, in inland waterways, and in the vicinity of river runoffs (Apel, 1987).Even in marginal seas and semienclosed basins swift flows may develop along semi-permanent circulation features (Robinson et al., 1999).However, the advances of operational oceanography show a great variability of the water flow at a wealth of spatial and temporal scales.This is indicated by both ocean drifter data -which are however affected also by wind (Maximenko et al., 2012), satellite altimetryfor just the geostrophic component of the currents (Pascual et al., 2006), and model computations -which capacity to represent mesoscale variability depends among others on spatial discretisation (Fu and Smith, 1996;Sandery and Sakov, 2017).More recently, even animal-borne measurements are employed to the end of characterising ocean currents -especially in the polar regions (Roquet et al., 2013).In view of applications, capturing such a complexity is mandatory for contributing to the value chain of ocean data (She et al., 2016).
The impact of ocean currents on navigation can be addressed from several viewpoints.Analytical changes of vessel heading for achieving least-time routes through a dynamic flow were derived from first principles by Techy (2011).However, they just hold for point-symmetric flow fields and constant vehicle speed with respect to the flow, which is definitely not the case for ship speed loss in waves.
An exact method based on the level set equation was developed by Lolla et al. (2014) and it is able to deal with generic dynamic flows and not constant vehicle speeds through the flow.It is based on two-step differential equations governing the propagation of the reachability front (a Hamilton-Jacobi level-set equation) and the time-optimal trajectory (a particle backtracking ordinary differential equation).The level set approach was extended to deal with energy minimisation by Subramani and Lermusiaux (2016) showing the potential of intentional speed reduction in a dynamic flow.This method appears to be quite promising, though it has not yet been embedded into an operational service.
Other mathematical techniques were reviewed in the introduction of Mannarini et al. (2016a) and some will be mentioned in this manuscript's Sect.3.1 for the sake of verification of the new numerical results.
In the latest edition of the World Meteorological Organization's guide to marine meteorological services, ocean and tidal currents are considered to be a key variable in the management of vessel fuel consumption (WMO-Secretariat, 2017).
The International Maritime Organization (IMO) recommends to avoid "rough seas and head currents" among the ten measures within the Ship Energy Efficiency Management Plan, or SEEMP (Buhaug et al., 2009).The SEEMP has come into force since January 2013, applies to all new ships of 400 gross tonnes and above, and is one of the main instruments for the mitigation of the contribution of maritime transportation to climate change (Bazari and Longva, 2011).

New contribution
The above recognition of literature shows that the question of the impact of sea or ocean currents on navigation, despite its classical appearance, is still open.In fact the available results are hardly comparable because: i) they are not validated versus exact solutions; ii) with some exception, they do not declare the computational performance; iii) generally, their model source codes are not openly accessible; iv) they are limited to case study analyses on a specific date, without an assessment of seasonal and geographical variability of quantitative conclusions; v) they generally cannot account for both waves and ocean currents.
All these considerations motivated a development of VISIR ship routing model as documented by the present paper, which is organised into three major sections: The theoretical framework for inclusion of currents into the model is presented in Sect.2; the verification of the numerics and computational performance is shown in Sect.3; the case-studies, comprehensive of an assessment of seasonal and geographical variability, are provided in Sect. 4. Finally, the concluding remarks in Sect. 5 precede the statement of the policy of availability of the model source code and input datasets.In App.A the main incremental changes of VISIR-I.b are documented.
We use throughout this manuscript the words "track" or "trajectory" for indicating a set of waypoints joining two given endpoints or harbours, in relation to departure on a given date, and the words "route" or "crossing" when there is no reference to a specific departure date.Furthermore, the shortcuts "w" for computations accounting for just waves and "cw" for both ocean currents and waves are used in the following.results that do not refer to any specific ship feature.The term "flow velocity" is used for referring to the velocity resulting from either ocean surface current, tidal current, and nonlinear mass transport in waves (Stoke's shift), or their composition.Also, when not otherwise specified, the qualification "over ground" is assumed for both speeds and courses.

Linear superposition
Assuming that a linear superposition principle holds for vehicle and horizontal flow velocity, the vector Speed Over Ground (SOG) is given by where F is the vehicle Speed Through Water (STW) and w the flow velocity.The symbol F is a reminder of the fact that such speed, due to energy loss in waves, in general is a function of both vehicle propulsion parameters and ocean state, cf.Mannarini et al. (2016a, Eq.21).
Eq. 1 is a "no-slippage" condition: the vehicle is advected with the flow.The rationale for this assumption is the experimental observation that ocean drifters (including vessels) adjust very quickly -i.e. in less than one minute -their speed to the flow (Breivik and Allen, 2008).At the present level of approximation, such adjustment is instantaneous (as no second derivatives of x appear in Eq. 1) and independent of vessel displacement (no vehicle mass in Eq. 1).Also Bijlsma (2010) and Techy (2011) Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.
in their optimal control methods and Zamuda and Sosa (2014), as a kinematic basis of an evolutionary approach for describing gliders' motion, make the assumption of linear superposition of speeds.Fossen (2012, Eq.26), in a context of vessel motion control, defines STW or relative speed through linear composition of SOG and current velocity.
However, we note that the linear superposition principle in the form of Eq. 1 just refers to a surface flow and cannot accommodate a depth-dependent (horizontal) flow speed w(z).In that case, vessel speed relative to water should be computed from the balance between the overall drag by the fluid (Newman, 1977) and the thrust provided by the propulsion system.This fact could be relevant for large draught vessels, especially those sailing in stratified waters (where the vertical profile of water velocity may exhibit both magnitude and direction changes, cf.Apel (1987)).
Finally, aerodynamic drag on vessel superstructure is also neglected in Eq. 1.

Course assignment
Along vessel path, course over ground (COG) may need to be constrained for navigational reasons (traffic constraints, fairways, shallow waters, or any other reason for preferring a specific passage).Also, in the computation of an optimal path, the algorithm (such as a graph-search method) may resort to spatial and directional discretisation, which again is a form of course assignment.
Making reference to Fig. 1, if COG has to be along ê, then the vehicle vector velocity must satisfy: where ô is a normal versor of ê.
In order to keep the course constrained as per Eq. 2, it is assumed that the shipmaster can act on the rudder(s) for modifying heading ĥ until COG satisfies Eq. 2 and then report the rudder(s) to midship.

Resulting kinematics
After defining the vector components of the water flow and making reference to Fig. 1, the flow projections along (ê) and across (ô) vehicle course (in either polar or rectangular representation) respectively are: where for both course ψ e and flow direction ψ w the nautical/oceanographic convention (i.e., "where-to" direction, clockwise from due North) is employed.Furthermore, the choice of orientation of the ô axis in Fig. 1 implies that a current bears to port whenever w ⊥ > 0. Linear superposition Eq. 1 and the course assignment condition Eq. 2 result into two scalar equations that, upon definition of an angle of attack δ of the ship's hull through the water (cf.Richardson (1997)): as the difference between the angle of vehicle heading (ψ s or HDG) and the COG, read with the unknown S g recognised as the vehicle SOG.Remarkably, Eq. 6a-6b could also be employed for determining ocean current vector w, given SOG, STW, course and heading.
As long as F is non null, δ is given by In presence of waves, F is reduced due to the wave-added resistance and can be obtained from a thrust-balance equation as in Mannarini et al. (2016a, Eq.14).Since F is always nonnegative, Eq. 7 implies that sgn(δ) = sgn(w ⊥ ).In particular, in case of a cross flow w ⊥ bearing to port, a clockwise change of vehicle heading is needed for keeping course, as in the example shown in Fig. 1.
Replacing δ into Eq.6a, SOG is obtained: Eq. 8 shows that the cross flow w ⊥ always (i.e., independently of its orientation) reduces SOG, as part of vehicle momentum has to be spent for compensating the drift.The along edge flow w ("drag") instead may either increase or decrease SOG.
Notice that the "cross" and "along" specifications refer to vessel course, differing from vessel heading by the (usually small) amount given in Eq. 7. Also it should be noted that the condition is not guaranteed in case of a strong counter-current.In a directed graph (as the one used in VISIR), violation of Eq. 9 along a specific edge would imply that the edge is made unavailable for sailing along that direction.
An equation formally identical to Eq. 8 was retrieved by Cheung (2017) in the context of flight path prediction, with wind replacing the ocean currents and plane true airspeed replacing STW.
Furthermore, both Eq. 7 and Eq. 8 hold if and only if If this is not the case, vehicle speed cannot compensate ocean current drift.We note that Eq. 10 is satisfied even in case of a vehicle drifting along the streamlines of the flow field without any steering (F = w ⊥ = 0).In that case, Eq.
Finally, by taking the module of both sides of Eq. 1 and approximating the l.h.s. with its finite difference quotient (thus leading to a first order truncation error), the graph edge weight δt is computed as where δx is the edge length.The weights δt are then employed for the computation of a time-dependent shortest path, through the same graph search method described in Mannarini et al. (2016a) and updated in this manuscript in Sect.2.4.

Navigationally safe graph
Due to non-convexity of the shoreline and presence of islands, the maritime space domain is not simply connected.This implies that not all graph edges correspond to navigable courses.Those unnavigable ones should therefore not be considered for the computation of the optimal paths.
In order to keep this into account, a graph pruning methodology had already been devised for VISIR-I.a(Mannarini et al., 2016a).Given the Under Keel Clearance (UKC) as the difference between sea depth and vessel draught, the method was based on checking both UKC > 0 and the absence of intersection between graph edges and coastline segments.In fact, just checking for a positive UKC is not enough.UKC is computed on a graph edge as the minimum value at the two edge nodes.However, UKC can still be positive for an edge crossing the shoreline with both nodes at locations with z > T .This required in VISIR-I.a to preliminary prune the shoreline-crossing edges, before the UKC > 0 condition was checked.
Since in a large ocean domain most of the edges do not intersect the coastline, a computationally more efficient procedure has been developed for VISIR-I.b.It consists in following three steps: i) Retrieve the indexes of edges within a small bounding box around each coastline segment; ii) Check edges within the bounding box for intersection with the coastline; iii) Create all edges in the selected domain, pruning just those -from ii) -which intersect the coastline.
The i) step can be performed in a constant time with respect to the size of the maritime domain by exploiting the fact that the graph is based on a structured grid.Furthermore, it can employ a lower-resolution version of the shoreline (cf.Sect.4.1.2) while the ii) step must employ a higher-resolution.
In the current version of the procedure, iii) step creates both land and sea edges.This is not an issue for navigation as, through the ii) step, such land edges are not connected to sea ones.However, the closest node to either of the track endpoints set by the VISIR user might be on land, making impossible to join it to the other endpoint.The user has then to manually move the endpoint.This is not a limitation for a not operational use, as in the present work, and will be fixed in future code versions.In VISIR-I.agraph nodes were just linked to all other nodes which can be reached via either one or two hops.In this work instead, a larger number of hops ν is allowed.This allows to increase the angular resolution ∆θ up to: The ν value is also called "order of connectivity" of the graph (Diestel, 2005).In Mannarini et al. (2018, submitted) the point is made that the numerical solution of the shortest path problem on a graph converges to the numerical truth as ν is increased in roughly inverse proportion to graph mesh spacing ∆ g 1 .
A limitation of the present graph generation procedure is the fact that multi-hop arcs (i.e.ν ≥ 1) intersecting landmass are being pruned, but those intersecting shoals are not.In order to improve on this, a further pruning step has been proposed in Mannarini et al. (2018, submitted).This effect is however more relevant in coastal and archipelagic applications than for the ocean crossings considered in the present work.
The computational cost of VISIR-I.b graph generation procedure is linear in the total number of edges (from step i) of the procedure above) within all the bounding boxes around the shoreline.For a given number of nodes, the computing time for preparing a graph of order ν then scales as O(ν 2 ) .

Time interpolation of edge weights
As in VISIR-I.a, also in VISIR-I.b edge weights are computed out of Eq. 12.
The shortest path algorithm is still derived from Dijkstra's one, which is a deterministic and exact method (Bertsekas, 1998).
The algorithm was made time-dependent under the assumption that no waiting times at the tail nodes are necessary, or "FIFO hypothesis" (Orda and Rom, 1990).Furthermore, a new option has been introduced in VISIR-I.b to perform time-interpolation of edge weights.In fact, the edge weights are no more kept constant between consecutive time-steps of the input geophysical fields (ocean currents and/or waves) but are estimated at the exact time the tail node is expanded by the shortest path algorithm.
In Mannarini et al. (2018, submitted) it was shown that the effect of time-interpolation can be relevant wherever the environmental fields rapidly change between successive time steps.This is likely the case for daily averages of the wave fields (Sect.4.1.4)which are employed for the case studies (Sect.4) of this manuscript.Orda and Rom (1990) stated that, under FIFO hypothesis, the worst-case estimate of the computational performance is, as for the static case, O(N 2 ), with the number N of graph grid points considered 2 .However, Foschini et al. (2014) made the point that, in presence of time dependent edge weights, the computational performance may degrade to become non-polynomial.The scaling of performance with time-interpolation is further investigated in Sect.3.2 through a few empirical test. 1 We here refer to a regular latitude/longitude mesh with ∆g spacing, distinguishing from its projection on planar coordinates, with a constant ∆y spacing and a ∆x depending on latitude.
2 The performance could be improved to O(N log N ) in a codification making use of binary heaps (Bertsekas, 1998).

Vessel modeling
VISIR-I.b vessel propulsion and seakeeping model is, with a minor update, the same of VISIR-I.a.It is shortly reviewed and updated in following Sect.2.5.1-2.5.2.Furthermore, under the hypothesis of constant Engine Order Telegraph (EOT), an estimate of the voyage energy efficiency is provided in Sect.2.5.3.

Vessel speed in a seaway
STW together with ocean current velocity determines SOG (Eq.1).SOG in turn determines the edge weights in the graph representation of the kinematical problem (Eq.12).STW depends on the vessel propulsion system (MANDieselTurbo, 2011) as well as on the energy dissipated through hydrodynamic viscous forces, aerodynamic forces, ocean gravity waves and waves due to water displacement (Richardson, 1997).It is however beyond the scope of this manuscript to develop a vessel propulsion and sea-keeping model more realistic than in VISIR-I.a(Mannarini et al., 2016a).
That model considered balance of thrust and resistance at the propeller, neglecting the propeller torque equation (Triantafyllou and Hover, 2003).In the resistance, a term related to calm water is distinguished from a wave-added resistance.The calm water term depends on a dimensionless drag coefficient C T that, within VISIR, is supposed to have a power-law dependence on vessel speed through water: C T = γ q (STW) q .For the wave added resistance, its directional and spectral dependence is neglected, and just the peak value of the radiation part is considered.The latter was obtained by Alexandersson (2009) as a function of vessel's principal particulars, starting from a statistical reanalysis of simulations based on Gerritsma and Beukelman (1972)'s method.Considering just radiation and neglecting the diffraction term, wave added resistance might be underestimated for vessels long with respect to the wavelength.

Vessel intact stability
In line with a IMO guidance (IMO, 2007), VISIR employs sea-state information also for performing a few checks of vessel intact stability.In Mannarini et al. (2016a) an ongoing research activity on this topic was noted.Specifically, at that time the development of "second generation" stability criteria had already been proposed by Belenky et al. (2011).A recent Terms of Reference for updating the IMO stability Code (IMO, 2008) has been published by the IMO Maritime Safety Committee (IMO, 2018c).
At present, VISIR includes checks of intact stability related to: parametric roll, pure loss of stability, and surfriding/ broaching-to at an intermediate level between IMO (2007) and the second generation criteria.Either intentional speed reduction (EOT<1, Tab. 1) or course change can be exploited by VISIR for fulfilling the stability checks (Mannarini et al., 2016a).Edges which, for a given EOT, violate stability are pruned before the shortest path algorithm is run.This way, it is ensured that the optimal track preserves vessel intact stability.
In this respect, in VISIR-I.b the sole update regards the actual values of the vessel parameters and the parametric roll stability check.
The new vessel parameters are suited for modelling a container ship and are listed in Tab. 4. These values result in a STW dependance on significant wave height as in Fig. 4a and resistances as in Fig. 4b.
For the parametric roll, the wave steepness criterion is generalised for vessels of L WL > 100 m by implementing the piecewise linear function of L WL given by Belenky et al. (2011, Eq.2.37).This means that Mannarini et al. (2016a, Eq.32) is replaced by where the critical ratio Σ is given by Since the stability changes are maximized for a ship length close to wavelength (Belenky et al., 2011, Sect.2.3.3), the Σ ratio represents also a critical wave steepness.Thus, Eq. 15 implies that it reduces at larger wavelengths, making the check on loss of stability in rough seas more severe than within the previous (VISIR-I.a)formulation.

Voyage energy efficiency
In this subsection the impact of track optimisation on voyage energy efficiency is estimated.The third IMO GHG study estimated the share of emissions from international shipping in 2012 to be some 2.2% of total anthropogenic CO2 emissions (Smith et al., 2014).According to the EDGAR database, emissions from international shipping in 2015 were larger than the quota of two Countries such as Italy and Spain altogether (JRC and PBL, 2016).
In line with the United Nations Sustainable Development Goal 133 , an initial GHG reduction strategy was approved by IMO in April 2018 (IMO, 2018b).It is layered into three levels of ambition, with the second one being "to reduce CO2 emissions per transport work, as an average across international shipping, by at least 40% by 2030, pursuing efforts towards 70% by 2050, compared to 2008".Furthermore, an implementation through short-, mid-and long-term measures is envisaged.The short-term measures include development of suitable indicators of operational energy efficiency.
IMO had previously introduced the Energy Efficiency Operational Indicator or EEOI as the ratio of CO 2 emissions per unit of transport work (IMO, 2009b).Depending on vessel type, there are several possible definitions of transport work.We here restrict ourselves to a cargo vessel carrying solely containers, for which transport work is defined as deadweight (DWT) times where the C F is a conversion factor from fuel consumption to mass of CO 2 emitted, s is the specific fuel consumption, P is the engine brake power and T the sailing time.Variations of P are allowed by VISIR algorithm, Sect.2.5.2, while s is assumed to be a constant.
If a track is plied at a constant P (i.e., EOT=1), the emissions are then proportional to T and the EEOI ratio ρ β,α of two tracks between same endpoints and sailed with same DWT is given by where the subscripts label the β track being compared to the α track.Eq. 17 shows that ρ β,α is the inverse ratio of the average speeds along the β and α tracks.The EEOI relative change of β to α track is then given by If the average speed in the β track is higher than in the α track, then −∆(EEOI) β,α > 0, i.e. a EEOI saving is achieved.
Depending on the subscripts α and β, different types of −∆(EEOI) β,α will be computed in Sect.4.4.4 for analysing the benefit of the optimal tracks.A not constant EOT is accounted for by the computation.However, for the EOT=1 limiting case, following general properties can be established: i) If vessel stability checks (Sect.2.5.1)do not lead to any diversions, the mean speed along the optimal track is never lower than along the shortest distance track.Thus, related EEOI savings are always non negative, −∆(EEOI) β,g ≥ 0; ii) Since currents can be either advantageous or detrimental to SOG (Eq.8), savings of the optimal tracks of cw-type can have any sign with respect to optimal tracks of w-type, −∆(EEOI) cw,w 0.
3 Verification and Performance VISIR-I.b path kinematics described in Sect. 2 is employed for the numerical computation of optimal paths on graphs.In this section, an assessment of VISIR-I.b numerics is provided by means of verification vs. analytical benchmarks (Sect.3.1) and test of its computational performance (Sect.3.2).

Analytical benchmarks
For the verification, VISIR-I.b includes an option for being run employing in input -in place of fields from data assimilative geophysical models (which will be described in Sect.4.1) -synthetic fields, leading to analytically known least-time trajectories or "brachistochrones".The rest of the processing (generation of the graph, evaluation of the edge weights, computation of the shortest path) is identical for both synthetic and modelistic environmental fields.However, as provided in Sect.3.1.1and Sect.3.1.2below, the synthetic fields are described in terms of linear coordinates.Thus, the spherical coordinates of the graph nodes are first linearised via an equi-rectangular projection.

Waves
The least-time route in presence of waves is computed by VISIR assuming that waves affect the speed through water of the vessel, Sect.2.5.1.For a static wave field, this leads to a STW not explicitly dependent on time.Thus, the least-time trajectory problem can be formulated in terms of a variational problem.
Analytical solutions are available for a subclass of these problems where STW depends on just one of the spatial coordinates (Morin, 2007).In particular, if speed through water F depends on the square root of position as in and the trajectory initial point is at y = 2R, the least-time trajectory is given by (an arc of) cycloid with R and g parameters determining length and acceleration, respectively (Broer, 2014;Jameson and Vassberg, 2000).The cycloid presents a cuspid at the initial point as, along a brachistochrone, the region with F = 0 has to be quit the sooner.The rest of the trajectory corresponds to refraction within layers of increasing speed or decreasing wave height, according to Snell's law.
The cycloidal benchmark was already exploited in Mannarini et al. (2016a).Thereto, the numerical error of VISIR-I.a in trajectory shape and duration was ascribed to the limited angular resolution (a graph with ν = 2 was employed).
For VISIR-I.b, we compute graphs of higher connectivity (Sect.2.3), allowing to approach the cycloidal benchmark more closely.The results are provided in Fig. 2a and Tab. 2. A relative error of less than 1 per mil in T * can be reached acting on just graph connectivity.This improves by about one order of magnitude on the accuracy of VISIR-I.a.
The cycloidal solution exploits the fact that a functional of the spatial coordinate is minimised under some necessary conditions provided by the Euler-Lagrange equations (Vratanar and Saje, 1998).The hypotheses leading to these equations are not satisfied for the more general case where the integrand of the functional explicitly depends on time.Instead, an assessment of VISIR solution in time-dependent waves was performed by comparison to the numerical results of an exact method based on partial differential equations (Mannarini et al., 2018, submitted).However, verification of VISIR with time-dependent fields versus an analytical benchmark is possible in the absence of waves and presence of currents, as described in following Sect.3.1.2.

Currents
This optimal control formalism provides the framework for computing extremals of a functional depending explicitly not just on spatial coordinates but also on time (Pontryagin et al., 1962;Bijlsma, 1975;Luenberger, 1979).It is based on the fact that the trajectory is controlled by a group of variables, for which an additional relation ("adjoint equation") holds.A variant of this approach, the Bolza problem, was employed for computation of optimal transatlantic tracks with a time-dependent STW by Bijlsma (1975).Due to topological issues, there are unreachable regions of the ocean, and the method involves guessing the initial vessel course, which may hinder the implementation in an automated system.Another variant is the approach by Perakis and Papadakis (1989), which accounts for a delayed departure time and for passage through an intermediate location.
However, its outcome is limited to finding just spatially local optimality conditions.
Instead, several benchmark trajectories are provided by Techy (2011) based on Pontryagin's minimum principle (Luenberger, 1979) and employing vehicle heading as a control variable.In particular, in presence of currents, and for a constant speed F relative to the flow (analogous of STW in the nautical case), an analytical relation between vehicle heading (which is the control variable) and vorticity of any (point-symmetric) current field is demonstrated.The field is given by: where both Γ and Ω may depend on time.For the case study (Techy, 2011, Example 3), the start and end points are set at the side of one equilateral triangle, which third vertex is at the flow origin (x = y = 0).Finally, the duration T * of the least-time trajectory is retrieved through an iteration on initial heading.
Fig. 2b displays VISIR.bsolution of problem Eq. 20 for a case where Γ is a non null constant (divergent flow) and Ω (one half of the vertical vorticity) linearly changes in time as per parameters of Tab. 2. Resulting optimal trajectory changes its curvature, swinging on both sides of the geodetic track, which is crossed at about one third of its length, cf.Techy (2011, Fig.12).The elongation of the swinging is quite small, with the optimal trajectory differing from the geodetic by less than 1% in length.This poses a challenge to the numerical solver on the graph, as many and accurate course variations are required on a short distance.Thus, it is not surprising to find that graph mesh spacing ∆ g is critical for achieving convergence, even more than the graph order of connectivity ν.However, this only holds if time-interpolation of edge weights (Sect.2.4) is employed.
Otherwise, no significant improvement can be achieved, as T * errors in Tab. 2 demonstrate.With VISIR-I.b, a minimum error of about 1.3% in T * is obtained for the graphs employed.

Computational performance
The computational performance of the new model version VISIR-I.b is here assessed.
The major changes in the source code with respect to the already published version (Mannarini et al., 2016a), are summarised in Tab.A1.All the computations for collecting the performance data were run on an iMac (Processor: 3.5 GHz Intel Core i7; RAM: 32 GB 1600 MHz DDR3).Results are displayed in Fig. 3. There, the number of degrees of freedom (DOF) is given by the product N t A of the number N t of time steps and the number A of graph edges.A in turn depends on the number of grid points N comprised within the geographical region selected and on the order ν of the graph.Jobs of DOF varying over more than four decades are considered, corresponding to various combinations of graphs order ν and mesh spacing ∆ g .Without time-interpolation, the optimal path algorithm is about a 8 times faster, Fig. 3c.Furthermore, in the same panel the computational overhead from use of ocean currents on top of waves is assessed.There is no overhead for the shortest path computations, as they employ in input a set of edge weights of the same size for both cases.Instead, edge weight values are determined through the specific environmental fields employed (waves alone or also currents).Thus, the preparation of the denominator in Eq. 12 causes an overhead for the total job, which is up to 30% for the sampled DOF range.Starting from 3×10 8 DOF (second larger job in size), a rise in the overhead is observed.In order to understand its origin, the computer RAM allocation is monitored in Fig. 3b-d.
The peak of the RAM is reported in Fig. 3b.Below 2 × 10 6 DOF, peak RAM is dominated by the reading of the input environmental fields and does not depend on the actual graph employed for the shortest path computations.Above that threshold, a rise of peak RAM is observed which, beginning at about 2 × 10 8 DOF, eventually saturates.Here, the limit of the computer physical memory is approached, which leads to swapping and to a degradation of performance, as already observed in Fig. 3c.
This is even more apparent in Fig. 3d, where the ratio of peak RAM for the cw-to w-type computations is displayed.Peak RAM allocation occurs -for large enough jobs -for edge weights preparation, prior to the run of the shortest path algorithm.
There is up to 50% extra RAM which needs to be allocated in case ocean currents are considered.In fact, the environmental scalar fields to be considered are five (significant wave height, direction, and peak period; zonal and meridional current), while the latter two are not employed for the w-type computations.Thus, at 2 × 10 8 DOF a sudden drop of the peak RAM ratio is recorded, as the allocation for the cw-case approaches the physical RAM while, for the w-case, it is still significantly lower than such a limit and can further grow.
Out of Fig. 3d it is possible to define a "computational efficiency region", for VISIR jobs with DOF lower than the one leading to the the drop observed in Fig. 3d.In fact, the computations in subsequent Sect. 4 are performed on a cluster with a RAM of 64 GB, allowing to operate in the efficiency region even for DOF as large as 5 × 10 8 .

Case studies
In this section, VISIR-I.b capacity to deal with both dynamic flows and sea state fields in realistic settings is demonstrated using ocean current and wave analysis fields from data-assimilative ocean models.
The section is organised into a presentation of the environmental fields employed for the computations, Sect.

Bathymetry
The GEBCO 2014 bathymetric database4 (Weatherall et al., 2015) is employed in VISIR-I.b.Its spatial resolution is 30 arcsec or 0.5 nmi in the meridional direction.

Shoreline
The Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG5 ) of NOAA (Wessel and Smith, 1996) is employed in VISIR-I.b.There are five versions (c, l, i, h, f) of the database, with a resolution of about two hundred meters in the best case.Depending on the geographic domain, VISIR-I.b employs different versions of the GSHHG for the generation the graph (Sect.2.3).This is for limiting the generation time in case of jagged coastline, such as in archipelagic domains. in VISIR-I.a for sailboats (Mannarini et al., 2015).Wind directly affects also motor vessels through an added aerodynamic resistance and a heeling moment, which are mostly important for vessels with a large superstructure, such as passenger ships (Fujiwara et al., 2006).This will be considered for future VISIR developments.
For the moment, a NOAA-Ocean Prediction Center review of marine weather6 is employed for describing the synoptic situation affecting the ocean state during the periods of the case study of Sect.4.3.Also, an archive of surface analysis maps7 is considered.

Waves
Wave analyses are obtained through CMEMS8 from the operational global ocean analysis and forecast system of Météo-France, based on third generation wave model MFWAM (Aouf and Lefevre, 2013).
It employs optimal interpolation of significant wave height from Jason 2 & 3, Saral and Cryosat-2 altimeters.The model also takes into account the effect of currents on waves (Komen et al., 1996;Clementi et al., 2017).To that end, surface currents The spatial resolution is 1/12 • (i.e. 5 nmi in the meridional direction).Three-hourly instantaneous fields of integrated wave parameters from the total spectrum (Spectral significant wave height, Mean wave direction, Wave period at spectral peak are averaged in a preprocessing stage based on "cdo dayavg" 9 into daily fields.Neither Stoke's drift nor the partitions (wind wave, primary swell wave and secondary swell wave) are employed yet.This is due to the computer RAM needed for storing and processing the fields for generating the edge weights (Sect.3.2).Due to a much larger fetch, the impact of swell is estimated to be more significant in the Southern than in the Northern Atlantic Ocean (Hinwood et al., 1982).
The wave dataset name is GLOBAL_ANALYSIS_FORECAST_WAV_001_027 10 ; the product validation is provided by a companion document 11 .The datasets have been downloaded from CMEMS at least 14 days after their date of validity, ensuring that the best analyses are employed.

Currents
Ocean currents are obtained through CMEMS from the operational Mercator global ocean analysis and forecast system, based on NEMO v3.1 ocean model, (Madec, 2008).
The dataset name is GLOBAL_ANALYSIS_FORECAST_PHY_001_024 12 ; the product validation is provided by a companion document 13 .The datasets have been downloaded from CMEMS (http://marine.copernicus.eu/)at least 14 days after their date of validity, ensuring that the best analyses are employed.

VISIR settings
For the results shown in this section, optimal tracks are computed on a graph with order of connectivity ν = 8 (cf.Sect.

Individual tracks
We first consider a transatlantic crossing in the northern Atlantic Ocean, between Norfolk, at the mouth of the Chesapeake Bay (37 • 02.5' N, 76 • 04.2' W) and Algeciras, just past Gibraltar Strait (36 • 07.6' N, 5 • 24.9' W).Both East-and Westbound tracks are considered, Fig. 5.
First of all, we note that the geodetic (or least distance) track is northwards bent, as it is to be expected from an arc of GC of the Northern hemisphere on an equi-rectangular projection.The Northern edge of the track is flattened due to the finite angular resolution of the graph: ∆θ ≈ 7.1 • from Eq. 13.Nevertheless, as Tab. 5 reports, the error made by VISIR in the length of the geodetic route is a few permil.This is comparable to the accuracy of the function for the computation of distances on the sphere (employed in VISIR) compared to the ellipsoidal datum (more accurate, but slower).

Meteo-marine conditions
The synoptic situation in the northern Atlantic during the week following June 21st, 2017 (departure date for the Eastbound track) was dominated by the Azores High blocking descent of subpolar Lows to the middle latitudes.This led to relatively calm ocean conditions (significant wave height H s < 5 m) for most of the region involved in the Norfolk-Algeciras crossing.
In the week following February 16th, 2017 (departure date for the Westbound track) a Low with storm-force winds formed near (41 o N, 52 o W) and then moved N, influencing wave direction on 19th and 20th.On February 22nd another storm with waves of H s > 8 m developed at (37 o N, 58 o W).
As the currents are concerned, we note that the Eastern edge of the crossing is N of Cape Hatteras and, thus, N of the GS branch called Florida Current (Tomczak and Godfrey, 1994).

Track spatial and dynamical features
The topological and kinematical features of the optimal tracks of the case study are discussed in this subsection.

Tracks topology
Four different solutions for the optimal tracks of the USNFK-ESALG route the are seen in Fig. 5 (red lines).
For the Eastbound voyage accounting for just waves (w-type, Fig. 5a) the optimal track is quite close to the geodetic one.This is due to the absence of waves of relevant height along the trajectory during the crossing (about eight days, cf.Tab. 5).
Discontinuities are seen between significant wave height fields at consecutive time steps (vertical stripes separated by dashed lines).This is enhanced by the daily averaging of the original 3-hourly fields, cf.Sect.4.1.4.
As the optimal track is computed for the same departure date and direction but accounting for ocean currents too (cw-type), the solution is significantly modified, Fig.On the Westbound voyage of w-type (Fig. 5c) the optimal track takes diversions both S and N of the geodetic track.This longer path can be sailed at an higher SOG than the geodetic one because it skips both the storm in the North-Eastern Atlantic at ∆t = 1 − 4 days since departure and the the storm developing at ∆t = 6 days, at the latitude of the arrival harbour.
The optimal track for the same departure date and direction but cw-type (Fig. 5d) leads to yet another solution with respect to the w-type track.In fact, it sails N of the geodetic all the time.The speed loss due to the encounter with the storm at ∆t = 2 − 3 days is balanced by the speed gains due to a meander of the North Atlantic current encountered at ∆t ≈ 4 days at 44 o N and by the benefit of sailing a bit further away of the rough sea than the corresponding w-type track at ∆t = 5 − 6 days.

Tracks kinematics
In order to get a deeper insight into the results, in Fig. 6 a few kinematical variables are extracted along both the optimal and geodetic tracks, for both cw-and w-cases.
Starting from the Eastbound route, Fig. 6a, the SOG of the cw optimal track greatly differs from corresponding geodetic track.Thanks to the GS in fact, SOG gains by up to more than 4 kn are experienced in the first half of the trajectory.During the final part of the navigation (∆t ≈ 6.5 days), a SOG > 22 kn peak appears shifted in both tracks.This is the signature of the Atlantic jet past Gibraltar, which is encountered about 5 hrs earlier along the optimal track (cf.below Fig. 6c).Instead, the SOG does not appreciably differ when w-type optimal and geodetic tracks are compared.This is consistent with the spatial pattern seen in Fig. 5a.
The geodetic Westbound track displays heavy oscillations in SOG with two deep local minima at ∆t ≈ 3; 6 days (Fig. 6b).
They correspond to the two storms NE and SW of the track mentioned earlier.The SOG differs from the one along the geodetic track just at ∆t ≈ 1.5 − 3 days along the optimal track of w-type, and this is due to its initial northbound diversion.Starting from ∆t = 4 days both optimal tracks significantly differ from the geodetic one, with the cw one proving to allow the larger SOG in the second part of the crossing.
In Fig. 6c-d the ocean flow component w along vessel course (Eq.4a) is displayed.This quantity, together with its normal counterpart w ⊥ , determines, through Eq. 8, the value of SOG.The difference between the optimal and the geodetic tracks is noticeable for both East-and Westbound navigation.In Fig. 6c it is seen that the algorithm manages to encounter a nearly always positive (i.e.along the course) w , which even exceeds 4 kn at the end of the first day.It is apparent that the same w oscillations are retrieved in the SOG linechart of Fig. 6a for ∆t < 3 days and at the ∆t ≈ 6.5 days peak.For Westbound navigation, w is mainly positive (apart from the initial impact of the Atlantic jet before Gibraltar is passed) along the optimal track and mainly negative along the geodetic, which sails against the GS.At ∆t = 4 days a NW-bound meander of the North-Atlantic current is encountered, with a positive drag of up to 1.5 kn.
Finally, the angle of attack δ needed for balancing the cross flow w ⊥ (Eq.5) is displayed in Fig. 6e-f.The track-average of δ is nearly zero, its maximum value is of the order of 10

Safety of navigation
The stability constraints of Sect.2.5.2 were checked for.However, some of them do not result in any graph edge pruning during the actual transatlantic crossing of the vessel under consideration (cf.parameters in Tab. 4).In fact, pure loss of stability is not realised due to the threshold condition on significant wave height of Mannarini et al. (2016a, Eq.36) not reached.
Surfriding/broaching-to is not activated due to the condition on Froude Number never larger than the critical one for the wave steepness encountered (Mannarini et al., 2016a, Eq.42-43).Employing the generalisation discussed in Sect.2.5.2, parametric roll may instead occur for the present vessel parameters and the North Atlantic wave climate.
Also, it is found that on this specific route and departure dates, the voluntary speed reduction (Sect.2.5.2) is not activated by the algorithm.This mean that the tracks are sailed at a constant P and that the CO 2 emissions are linearly proportional to sailing time T * (Sect.2.5.3).Instead, for other routes in the Atlantic, this is not always the case, cf.Tab. 7.
Furthermore, all time-dependent edge weights along the optimal tracks fulfil the FIFO hypothesis (Sect.2.4).

Track metrics
Two simple metrics for summarising the kinematics of a track are here proposed: the optimal track duration T * and corresponding length L (not starred symbol, as such length is not object of optimisation).For the geodetic tracks, optimisation is instead performed on length L * and, unless safety constraints play a role in the actual optimal track, corresponding duration T is higher than T * .
L is sensitive to the geometrical amount of the track diversions, while T * reflects their kinematical impact.Such key metrics are reported in detail for both the geodetic and optimal tracks of both the East-and Westbound crossings in Tab. 5.The data also allow to distinguish the quantitative role of waves and currents and the amount of the track duration gains.For instance, it is seen that both East-and Westbound tracks lead to time savings ∼ 3% with respect to the geodetic track.However, for the former such a saving is mainly due to exploitation of currents, while for the latter to waves.
Concerning time gains, it is important to specify whether they refer to the geodetic track (∆T g ) or to an optimal track computed in presence of waves only (∆T w ).At this place, we observe that both Lo and McCord (1995) and Chang et al. (2013), not employing waves, just consider ∆T g .Also, in their case, the model region chosen for their track optimisation nearly coincides with the domain where the western boundary current under consideration is at strongest.This is at a difference with the case study presented in this section, which entails also the Eastern part of the ocean, where the influence of the Western boundary current is less noticeable.Thus, ∆T g gains due to currents reported in Tab. 5 are lower than those literature results, though possibly more realistic, because referring to full transatlantic crossings.

Track seasonal variability
In this subsection we address the question to what extent the seasonal variability of ocean state and circulation affects the variability of the optimal track of a given transatlantic crossing.
In order to address it, VISIR-I.b computations are carried out for departure dates spanning the whole calendar year 2017.
In particular, for each month, departures on six dates (1st, 6th, 11th, 16th, 21st, and 26th) are considered, leading to 72 dates per year.This is meant to account for decorrelation of the ocean current fields after a Lagrangian eddy timescale of about five days (Lumpkin et al., 2002).As waves are mainly driven by winds, which velocity is one order of magnitude larger than ocean velocities, the timescale for decorrelation of ocean state is expected to be even shorter.

Spatial variability
A direct visualisation of the annual variability of the track topology is shown in Fig. 7.
Each panel displays a bundle of trajectories relative to 72 departure dates.The extent of the diversions makes clear that the case study of Sect.4.3 is not even extreme.Instead, for both East-and Westbound tracks, the Summer and Autumn tracks are closest to the GC track.This is due to the fact that, in the Northern Atlantic Ocean, wave heights tend to be smaller in those seasons and, consequently, both vessel speed losses and relative kinematic benefits from diversions, are smaller.Some tracks are found to sail quite inshore into the Canadian coast, see a related comment in Sect.4.5.4.
The general impact of ocean currents on Eastbound tracks is that the bundle of tracks squeezes and shifts S in the vicinity of the GS proper (W of 67 • W).On a few dates (mainly in Winter and Spring) this is not the case, as storm systems happen to cross the location of the GS.For the Westbound tracks, accounting also for currents just adds a small perturbation to the wave-only tracks, without dramatically changing their topology.
It should be stressed that the computed spatial variability heavily depends on how ship energy-loss in waves is parametrised, cf.Sect.2.5.1.In fact, wave-added resistance determines vessel STW for any given sea state and, thus, how profitable a diversion is to the end of avoiding a speed loss.

Evolution lines
While the paths of the tracks displayed in Fig. 7 convey the information about the spatial variability and its seasonal dependence, they fail to inform about vessel kinematics along the tracks.To this end, an alternative visualisation is proposed in Fig. 8.
There, following a practice used in track anomaly detection (Zor and Kittler, 2017), cumulative sailed distance is displayed vs. Furthermore, in presence of currents, the slope can exceed the one relative to navigation at SOG equal to maximum STW.This is due to speed superposition per Eq. 8 and is apparent for some of the Summer tracks in the panel relative to the Eastbound tracks, Fig. 8c.
Finally, the envelope of the evolution lines along the geodetic tracks is displayed as a grey etched area.This makes the kinematical benefit of the optimal tracks apparent, as the optimal tracks can be sailed at an higher SOG (coloured dots are generally left of the grey areas), resulting in a shorter duration of the voyages.

Scatter plots
In order to reduce and better analyse the information contents of Fig. 8, the compound metrics T * and L can be employed, which are reported in a Cartesian plane in Fig. 9.
Such a plane contains a strictly forbidden region, left of L = L GC , which is the length (on the graph) of the GC arc connecting route endpoints.The straight line through the origin, which slope is V −1 max , generates another relevant partitioning of the plane.In fact, the region upper (lower) of this line corresponds to tracks sailed at an average speed lower (higher) than V max .
We first focus on Eastbound tracks.The distribution for w-type tracks is given in Fig. 9a.As expected, they are all comprised within the region above the T * = L/V max line.This is due to involuntary speed loss in a seaway, which reduces the average speed to less than V max .As also currents are considered, Fig. 9c, the tracks can be faster and, for Eastbound navigation, some of them even attain the region where the average SOG is larger than V max .This generally occurs for Summer tracks, which experience a lower speed loss in waves.
For the Westbound tracks, Fig. 9b-d, the general picture differs by the following features: the region where the average vessel SOG is larger than V max is never attained and the distribution in the (L, T * ) plane roughly maintains its pattern among the wand cw-type results.
These findings are also mirrored in the Pearson's correlation coefficient R P between T * and L. While for the Westbound tracks R P is nearly unchanged (Fig. 9b-d), it decreases substantially between Fig. 9a-c.This is due to the fact that most Eastbound tracks, independently of their duration, require a significant diversion to exploit the GS proper.This in turn reduces the correlation between T * and L.
The dots relative to the tracks selected for the featured analysis of Sect.4.3 are circled in Fig. 9.For the Eastbound crossing, a transition into the efficiency region is seen when comparing the w-to the cw-tracks.

EEOI savings
For assessing the benefit of track optimisation in terms of voyage energy efficiency, in Fig. 10 the monthly and annual variability of the EEOI savings are displayed.
In reference to Sect.2.5.3, specific fuel consumption s is taken to be a constant, while engine brake power P is allowed to vary as EOT is selected by the optimal routing algorithm (cf.Sect.2.5.2).
With the notation of Eq. 18, EEOI savings of the tracks considering both ocean currents and waves (β =cw) are computed with respect to either the geodetic track (α =g, Fig. 10a-b) or the wave-optimal tracks (α =w, Fig. 10c-d).
For the Eastbound route, −∆(EEOI) cw,g exhibits a clear seasonal cycle, with a peak of the monthly-mean value in Winter.
However, Winter intra-monthly variability exceeds the amplitude of the seasonal cycle.For the Westbound route, these trends are still observed, but both the seasonal cycle and the intra-monthly variability are less regular.
Furthermore, in Fig. 10c-d it is seen that the monthly-mean value of −∆(EEOI) cw,w is larger for the Eastbound route, since it can benefit from advection by the GS.Peak values of −∆(EEOI) cw,w are found in Summer months, when ocean state is calmer and thus the relative contribution of currents is the prevalent one.
Thus, the magnitude and location of the GS is critical for voyage energy efficiency along this route in Summer.In this respect, Minobe et al. (2010) found from satellite altimetry data that the seasonal cycle of the geostrophic component of the GS is weak both in meridional position and near-surface velocity.Instead, the simulations by Kang et al. (2016) show a seasonal cycle of the mean kinetic energy of the GS proper, with a relative maximum during Summer.Berline et al. (2006) analysed the GS latitudinal position at 75-50 • W from model re-analyses, finding that inter-annual and seasonal variability dominates upstream and downstream of 65 • W respectively.

Ocean-wide statistics
For demonstrating the generality of the VISIR-I.bcode and for assessing the potential EEOI savings depending on various wave and ocean circulation patterns, optimal routes are computed in the whole Atlantic Ocean.This requires among others that graph, shoreline, bathymetry, and environmental datasets of waves and ocean currents are made available to VISIR for wide enough regions for accounting for the spatial variability of the tracks.
Transatlantic crossings may in some cases exceed ten days, ruling out the use of wave forecast model outputs, which are limited by the availability of related atmospheric forcing fields.In fact, the maximum lead time of CMEMS products is ten days for ocean current forecasts but just five days for wave forecasts (cf.Product User Manuals cited in Sect.4.1).To our knowledge, at the moment ECMWF 14 runs a global wave model based on WAM with ten days lead time, but at a lower spatial resolution (1/8 • ) and without an open access policy, while NCEP 15 runs a model based on WW3 on various grids and lead time of 7.5 days 16 .
The unavailability of long enough forecasts can be addressed by either re-routing or use of supplementary information.
Re-routing or re-planning is the dynamic update of the optimal track as new information (forecast) is made available (Stentz et al., 1995;Likhachev et al., 2005).Corresponding solution is sub-optimal, as the initial routing choices are unrecoverable and may compromise the reaching of a global-optimal solution.An example of use of supplementary information instead has been 14 https://www.ecmwf.int/en/forecasts/datasets/set-ii 15http://polar.ncep.noaa.gov/waves/hindcasts/prod-multi_1.php 16https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfsproposed by Aendekerk (2018).Thereto, a "blending" of climatologies and geometrical information is employed as a surrogate for missing forecasts at large lead times.
In a not operational mode, the unavailability of forecasts is not critical.In that case, analysis fields can be employed, enabling a better reconstruction of the environmental state.A product derived from analyses may be quite useful for scenario assessment, while the uncertainty associated to forecast (Bos, 2018) may complicate its usefulness to that end.Analysis fields of waves and ocean currents are employed throughout the present manuscript.
For 9 ordered couples of harbours from the list in Tab.6, 72 tracks relative to year 2017 are computed.Furthermore, both sailing directions and both w and cw cases are considered, leading to the computation of 288 tracks/route/year.This sums up to the computation of more than 2,500 tracks in the Atlantic Ocean with the same VISIR-I.bcode version.
With the help of Tab. 7 and Fig. 11, some general results can be summarised as follows: a) EEOI savings are dominated by waves in the Northern Atlantic, with a not negligible contribution from currents.At the Equator, currents are the main reason for EEOI saving.In the Southern Atlantic, the largest savings are computed, and they are mainly ascribed to waves; b) Routes mostly impacted by ocean currents exhibit a large reduction of the correlation coefficient R P when comparing w-to cw-type scatter plots of track duration vs. track length; e) The FIFO hypothesis is not satisfied in just a tiny number of edges not employed for the optimal tracks.This backs the use of a time-dependent Dijkstra's algorithm as in Sect.2.4; d) Intentional vessel speed reduction (EOT<1) occurs in just three routes and for a quite limited fraction of their track waypoints.This backs the approximation done in Sect.2.5.3 for estimating the relative EEOI savings; e) Maximum ROT is never larger than about 20 • /min, which should be feasible for vessel manoeuvring, given typical accuracy of 1 • /min for onboard ROT Indicators (IMO, 1983).ROT is computed in VISIR-I.b from HDG changes.HDG differs from COG due to cross currents, Eq. 5. Thus, large values of ROT are experienced in routes crossing strong or variable currents.
In the following 9 paragraphs, some route-specific results are discussed.In the Supplementary Material of this manuscript related figures are published.

Buenos Aires -Port Elizabeth
The geodetic track is bent southwards in Mercator projection.The (Northern Hemisphere) Winter tracks are closer to the geodetic, while Summer ones exhibit the larger diversions.This route is characterised by the highest impact of waves on energy efficiency savings.This can be ascribed to the strength of the Antarctic circumpolar winds, causing large waves in the Southern Ocean (Lu et al., 2017).The role of currents is instead about 1%, with a stronger contribution from the Benguela current for Eastbound crossings.This is generally due to the avoidance of the Agulhas current past Cape Town.Biscay, no matter if ocean currents are accounted for or not.This is due to the activation the parametric roll safety constraint (Sect.2.5.2), as the encounter period of waves is about half the natural roll period T R of the vessel, Tab. 4. This occurs just for the track leaving from Rotterdam, as waves are encountered at a lower frequency on the other sailing direction.

Miami -Panama
The spatial variability of this route is dominated by currents, as waves from sub-polar Lows are not relevant in the Caribic region.The bundle shows a waist W of Cuba (21 • 52' N, 85 • 00' W), a point through which all optimal tracks but one sail.In fact, on Sept.11th, 2017 the track leaving Miami is affected by large waves in the Gulf of Mexico generated by the transit of Irma hurricane17 .In that case, the sea state, together with a local intensification of the GS in the Florida straits, leads to an optimal track sailing E of Cuba.

Boston-Miami
This route is heavily influenced by the Florida current.The Northbound tracks tend to align with the ocean flow.The Southbound tracks (sailing against the main flow) split into two sub-bundles, W and E of the Florida current.The Western sub-bundle is populated by mainly Winter tracks.In fact, these tracks sail more inshore, avoiding rough ocean state and, thus, reducing speed loss in waves.

Conclusions
VISIR ship routing model and code have been updated to version I-b.Optimal tracks can now be computed in presence of both time-dependent ocean currents and waves.Vessel interaction with currents is described in terms of new equations which are validated by means of an analytical benchmark.Furthermore, in order to represent vessel courses at a higher accuracy, the previous model version has been improved with respect to the capacity of computing graphs of higher order of connectivity, keeping the shoreline into account.The computational cost of the new model version is also assessed, and inclusion of ocean currents leads to a total CPU time overhead not exceeding 30% for realistic computations (Fig. 3c).
While the code of VISIR-I.a was tested by its operational implementation in the Mediterranean Sea (Mannarini et al., 2016b), the robustness of VISIR-I.b has been proven through computation of more than 2,500 tracks via the same model code version, ranging nearly any subdomain of the Atlantic Ocean.
Several routes are considered, mapping the variability of the optimal tracks along a full calendar year (2017).Both spatial and kinematical variability of the tracks are accounted for, through various types of diagrams.Optimal exploitation of ocean currents may lead in some cases to average speeds greater than maximum vessel speed in calm water (cf.Fig. 8-9).Finally, a standard voyage efficiency indicator (EEOI, introduced by the International Maritime Organization) is employed for highlighting the contribution of ocean currents and waves to the efficiency of the voyages.In some cases, EEOI relative savings in excess of 5% (annual averages) and 10% (monthly averages) can be reached (Fig. [10][11].However, the intra-monthly, seasonal, and regional dependence of these results is quite strong, and this study provides one of the first attempts to quantify it.Also, it should be recalled that these results depend on the actual parametrisation of wave-added resistance, which is still formally the same of Mannarini et al. (2016a).These quantitative assessment of EEOI savings through path optimization may be considered in relation to the ongoing discussion at IMO-level on comparing the effectiveness of several proposed methods for vessel emission savings (IMO, 2018a).
Furthermore, the analysis of the track dataset is simplified by means of metrics such as the optimal track duration and length, their Pearson's correlation coefficient, and the maximum rate of turn of vessel heading.The correlation coefficient carries a signature of ocean currents, which tend to make optimal track duration and its length less correlated to each other.Furthermore, the approximation of a FIFO network (Sect.2.4) is monitored and found to be satisfied to a very wide extent (Tab.7).Vessel EOT is allowed to vary (Sect.2.5.2), and the computation of the EEOI savings do account for it.However, intentional speed reduction is found to be a rare choice of the optimization algorithm.
We think that the major computational limitation of VISIR-I.a,b is its requirement on computer RAM allocation (Sect.3.2).
In fact, the code still requires preparation of all the time-dependent graph edge weights, ahead of the shortest path computations.
This presently impacts the capacity to describe the environmental state surrounding the vessel.For instance, in this work we averaged three-hourly wave fields to daily averages (Sect.4.1.4),we neglected other wave spectrum components (such as swell), nor we did account for the Stokes's drift contribution to the flow advecting the vessel.
On the other hand, it should be remarked that a more realistic representation of the marine state should correspond to a more accurate description of the mechanical interaction between it and the vessel, especially with reference to speed loss in waves and wind (Tsujimoto et al., 2013).Also, presence of sea ice and ECA zones may affect the optimal tracks.While the former effect may decrease in significance in a warmer planet, the latter has the potential shape more and more coastal traffic, as the new IMO global cap on sulphur contents enters into force (IMO, 2016).Progress on the representation of some of these model components is planned for future VISIR versions and will pave the way for end-to-end model evaluation exercises with respect to actually sailed trajectories.Table 5. Route length L (or L * for geodetic tracks) and optimal duration T * (or T for geodetic tracks) for tracks in Fig. 5-6.∆Tg is the relative duration change with respect to the geodetics; ∆Tw with respect to w-type optimal tracks.On the WGS-84 geoid, the length of the arc of GC between the endpoints is 3332.60 nmi, i.e. the numerical solution overestimates it by 0.3%.In a still ocean (no currents nor waves) the numerical geodetic would be sailed in 158 :28 :28 hrs by a vessel with Vmax as in Tab. 4. The second header line of the −∆EEOI columns specifies the type of tracks as in Eq. 18. track Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.

Following
the Paris Agreement (UNFCCC, 2015), climate change of anthropic origin is attaining increased attention at international and regulatory level.The Intergovernmental Panel on Climate Change recently published a special report addressing the Green House Gases (GHG) emission reduction pathway for complying with a global 1.5 • warming above pre-industrial levels.It was noted that this would require rapid and far-reaching transitions in energy systems and transport infrastructure (IPCC, 2018).

Fig
Fig.3adisplays both the cost of computing just the optimal track via the shortest path algorithm and the total job cost, since submission to saving of the results (rendering excluded).It is distinguished between the cases without and with time interpolation of the edge weights (Sect.2.4).The CPU time for the optimal track scales nearly linearly with the DOF.Instead, 4.1; a documentation of the principal VISIR model settings employed, Sect.4.2; a description of the results on individual tracks of given departure date, Sect.4.3, the analysis of their seasonal variability within a calendar year, Sect.4.4, and the extension of such analysis to several routes in the Atlantic Ocean, Sect.4.5.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.4.1 Environmental fields VISIR-I.b employs both static and dynamic environmental fields from official European and US providers.The static environmental datasets are bathymetry and shoreline.The dynamic datasets are waves and ocean currents.The specific fields employed are described in the following subsections.
Meteorological fields are not yet employed for computing VISIR-I.b tracks.So far, surface wind fields have just been used Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.from corresponding CMEMS product (see Sect. 4.1.5)are employed and used to force daily the wave model.The currents modulate wave energy and also cause a refraction of the waves propagation.
2.3)and mesh spacing ∆ g = 1/8 • .This graph resolution parameters are chosen to compromise between track accuracy (i.e.spatial and angular resolution) and computational cost of the numerical jobs.The computations refer to a container ship, which parameters are reported in Tab. 4 and which performance in waves (computed by the same method ofMannarini et al. (2016a)) is summarised in Fig.4.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.
5b.A diversion S of the geodetic track is computed by VISIR-I.b.It is instrumental in exploiting advection by the GS through velocity composition (Eq.8).Despite longer in terms of sailed miles, this track is faster than the geodetic one, Tab. 5. A closer look at Fig. 5b reveals that the optimal track averages between the locations of Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.opposite meanders of the first six oscillations of the GS proper, at 72-63 • W. Subsequent meanders, being prone to extrude filaments (and thus more stretched in the meridional direction), are followed less and less closer by the optimal track.
Fig. 6e-f are compliant with vessel manoeuvrability.The maximum module of the Rate Of Turn (ROT) of HDG is found to be 2.9 • /min for the East-and 1.5 • /min for the Westbound track of cw-type.These values are comparable to the IMO prescribed accuracy of 1.0 • /min for onboard ROT Indicators (IMO, 1983).Thus, heading fluctuations computed by VISIR-I.b for this route are feasible with respect to manoeuvrability.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-292Manuscript under review for journal Geosci.Model Dev. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.time elapsed since departure.Thus, the slower parts of each trajectory result in a smaller slope for corresponding segments of the track "evolution line".It can be seen that such slow segments are more frequent in Winter months and in the middle of the crossing, especially for Westbound tracks.This is due to larger speed losses in waves.

Figure 1 .Figure 2 .
Figure 1.VISIR-I.b directional conventions on top of the compass protractor.Vessel speed through water F, flow speed w, vessel COG ψe, vessel heading ψs, flow direction ψw, angle of attack through water δ = ψs − ψe.The longer (shorter) cathetus of the blue (red) triangle is the cross current magnitude w ⊥ = F sin δ, cf.Eq. 7. The displayed configuration refers to a vessel course assignment ψe = 25 • and implies a positive angle of attack δ = 21 • for balancing the drift due to a port-bearing flow w.

Figure 3 Figure 4 .Figure 6 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 3. a) CPU time for the total VISIR job (blue markers) and for just the computation of the time-dependent shortest path (red markers).The cw case only is considered.Dashed lines are fits of the model in Tab. 3. b) Peak RAM allocation during the jobs of a) panel, with a reference line at the total installed RAM.c) Ratio of CPU times between the cw and w cases and (just for optimal path) with and without time-interpolation.d) Ratio of peak RAM allocation of the cw to w type jobs.For panels a,b,d) both cases with (filled) and without (empty markers) time-interpolation.The DOF (Sect.3.2) of the time-dependent shortest path problems is displayed on the horizontal axis.
is the length of the geodetic track on the graph.∆ is a shortcut for the EEOI saving.The • operator denotes the annual mean, the • the mean annual value of the standard deviation.Corresponding values are given in %.The second header line specifies the type of tracks.The other columns contain: the number of tracks NE with intentional speed reduction and the maximum fraction of track waypoints affected (WP ) -for the w-type this figure is always 0 but for the ZAPLZ-ARBUE route, where it reads 1(0.4); the maximum rate of turn ROTM ( • /min); the number of non-FIFO edges F (neither of them is along the optimal track); the Pearson's coefficient RP between T * and L. The DOF varies from more than 5.4•10 8 of the ARBUE-ZAPLZ to about 2.5•10 7 of the USBOS-USMIA.

Table 1 .
Some nautical abbreviations employed in this manuscript.

Table 2 .
Summary parameters of benchmark case studies, cf.Fig. 2. Length scale L0 set by the track endpoint distance and time scale T0 = L0/Vmax are employed throughout (Vmax = 21.1 kn).Values in italics correspond to runs without time-interpolation of the edge weights, cf.Sect.2.4.Values in the last row of each group refer to the analytic solutions.

Table 3 .
Fit parameters for the data displayed in Fig.3a.The fit model is a • x b + c.For the optimal path data, c parameter is not fitted.

Table 4 .
Database of vessel propulsion parameters and principal particulars used in this work.The values of ∆ (ballast -scantling) and the maximum cargo capacity (2 500 TEUs) are not used in the computations and are provided just for the sake of reference.

Table 6 .
direction track type forcings L (or L * ) T * (or T ) Database of harbours.Coordinates refer to the graph grid point selected by VISIR.Wherever available, GRT is the annual throughput in year 2016 from Lloyd's (2018) and is used for sorting the entries.The other harbours are sorted alphabetically by International Seaport

Table 7 .
Database of routes.Lg

Table A1 .
List of main code changes of VISIR-I.b with respect to VISIR-I.a version, with indication of their use within this manuscript.