Ocean Engineering 195 (2020) 106706

Contents lists available at ScienceDirect

Ocean Engineering
journal homepage: www.elsevier.com/locate/oceaneng

Development of a numerical marine weather routing system for coastal andmarginal seas using regional oceanic and atmospheric simulations
Kenta Kurosawa a, b, *, Yusuke Uchiyama c, d, Taichi Kosako da

RIKEN Center for Computational Science, Kobe, Japan
Department of Atmospheric and Oceanic Science, University of Maryland, College Park, MD, USAc
Department of Civil Engineering, Kobe University, Kobe, Japand
Coastal and Estuarine Environment Department, Port and Airport Research Institute, Yokosuka, Japanb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Weather routing
A-star algorithm
ROMS
Kuroshio
Optimal route
Numerical simulations

The primary objective of this study is to develop a new weather routing system for vessel navigation in coastal
and marginal seas based on the A-star algorithm. The cost function inherent in the original A-star algorithm ismodified to account for oceanic and atmospheric conditions around the vessel of interest. Three options are
introduced to search for optimal paths with the shortest travel distance, shortest travel time, and minimal fuelconsumption. An avoidance algorithm for unsafe conditions is further incorporated to exempt any arbitrary areafrom the navigation solution. Furthermore, a compact ocean circulation model based on the Regional Ocean
Modeling System executable on typical laptop PCs is configured for a vessel-borne weather routing system; thesystem is successfully applied to evaluate the optimal vessel paths in the Seto Inland Sea, where high-frequency
tidal currents modified by complex topography make it essential to alter the vessel’s speed.

1. Introduction

factors such as weather and oceanic conditions, characteristics of ships,
and areas of territorial waters (Sen and Padhy, 2010; Tsou and Cheng,2013). WR differs from automotive navigation systems in many aspects.
For example, while automotive navigation systems usually consider only
traffic congestion and negotiable roads, WR must account for natural
phenomena such as oceanic currents, wind speeds, and wave heights.
Moreover, in the open ocean, there are very few obstacles except islands
and territorial waters in contrast to the land. Hence, because WR must
deal with an infinite number of routes consisting of vast numbers of links
and nodes but with fewer guideposts to find the optimal route, its
computational cost becomes drastically more expensive than that for
automobile navigation. Furthermore, it is necessary to design WR al­
gorithms properly not only for the shortest route (Efentakis et al., 2011)
but also for the safest route (Kosmas and Vlachos, 2012), the minimum
fuel route (Jung and Rhyu, 1999; Papanikolaou et al., 2016), and the
combinations of these factors (Szlapczynska, 2007; Andersson, 2015).
Therefore, previous studies have attempted to develop path-finding al­
gorithms while minimizing the computational cost of WR.
One of the more popular methods for finding an optimal route for WR
is the isochrone method proposed by Journ�
ee and Meijers (1980).
Although this method has been used for decades because of its
simplicity, the isochrones have to be manually determined a priori.

Reduction of carbon dioxide (CO2) emission has been a long-standing
issue for mediating global warming. Among others, fuel efficiency inship navigation is recognized as a fundamental factor that must be
remedied. Buhaug et al. (2009) presents that the proportion of CO2 fromship emissions in 2007 was estimated to be 1,046 million tons, ac­counting for 3.3% of the global CO2 emission during that year. For
economical operation in marine logistics, optimal routing to determine
on-time or shortest cruising path is indispensable. In particular, finding
the shortest path is often more essential to immediate excursions forfishing, rescues, and leaked contaminants, and to provide urgent evac­
uation from coastal disasters such as tsunamis and storms. However,vessel navigation is affected significantly by ambient oceanic and at­mospheric conditions that restrict the optimal paths. For instance,
opposing oceanic currents reduce the vessel’s relative speed or increase
fuel consumption, whichever the operator cares for more seriously(Trodden et al., 2015; Bouman et al., 2017). Therefore, weather routing(hereinafter denoted as WR), in which surrounding oceanic and atmo­
spheric conditions are taken into account in navigation, is highly
desirable for desirable vessel navigation.
In general, WR determines the optimal route by considering many

* Corresponding author. RIKEN Center for Computational Science (R-CCS), 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo, 650-0047, Japan.
E-mail address: kkurosaw@umd.edu (K. Kurosawa).
https://doi.org/10.1016/j.oceaneng.2019.106706
Received 3 July 2019; Received in revised form 30 September 2019; Accepted 7 November 2019
Available online 19 November 20190029-8018/© 2019 Elsevier Ltd. All rights reserved.

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

such algorithms directly for vessels sailing on coastal and marginal seas
(Smierzchalski and Michalewicz, 2000; Chen et al., 2013). The primaryreason is complexity of weather and oceanic conditions as well as
coastline topography, all of which are necessary for coastal WR. Thus,
there is room for improvement on WR in such areas. One of the moresignificant issues arises from a poor resolution of the oceanic and at­mospheric information released from agencies concerned. In general,these products have quite large spatial resolutions and large update
intervals. For example, a regional oceanic reanalysis/forecast product,
the Japan Coastal Ocean Prediction Experiment (JCOPE2; Miyazawa
et al., 2009) of the Japan Agency for Marine-earth Science and Tech­
nology (JAMSTEC), provides daily-averaged data at a lateral grid
spacing of about 10 km. This resolution is not sufficient in crowded
harbors and bays, where detailed navigation is crucial, because
high-frequency variability such as tidal currents characterize theambient oceanic condition most prominently. Moreover, data retrieval
prior to departure and subsequent onboard predictions are occasionallyindispensable because ample data retrieval through communications
satellites during the voyage may not be feasible. Therefore, this study
also configures a compact and computationally efficient regional ocean
model executable on typical laptop PCs aiming at a vessel-borne weather
routing system for suboptimal availability of weather and oceanic
reanalysis/prediction. The model is coupled with the WR system basedon the A-star algorithm to evaluate the optimal vessel paths. Our test bed
is the Seto Inland Sea (hereafter SIS), Japan, filled with an enormous
number of commercial carriers and small fishery and pleasure boats. In
the SIS, the high-frequency tidal currents modified by complex topog­
raphy are substantial to alter the vessel’s speed for optimal navigation.
Therefore, this study aims to achieve the following.

Hagiwara (1989) improved the isochrone method to be suitable forcomputers by implementing an automated isochrone detection algo­rithm. However, the isochrone method occasionally generates so-called“isochrone loops” making path finding impossible, and failing to avoid
crossing lands. Hence, De Wit (1990) and Calvert et al. (1991) intro­
duced dynamic programming, which uses a grid system that divides a
possible sailing region into several cells where each crossing point of acell boundary is used as a waypoint candidate. Some studies have madeefforts to conduct a series of detailed investigations to reduce the nav­igators’ workload by using a dynamic programming method for grid size
setting and design of the WR system (Wei and Zhou, 2012; Kim and Lee,2018). Although these methods have been developed elaborately to
account for oceanic and atmospheric conditions as much as possible,
their main drawbacks are computational inefficiency caused by a large
number of nodes, a requirement for large storage space, inflexibility, andslow computational speed.
With a gridded sea area, Dijkstra’s algorithm (Dijkstra, 1959) ingraph theory is well known to be practical for ship navigation. The al­
gorithm begins at a specific node and extends outwards to calculate theshortest path between any two vertices in a weighted graph (Fig. 1a).
Because the algorithm requires searches for all the nodes between theorigin and destination, it is computationally very expensive if the sea
area of interest is gridded at fine lateral spacing. By contrast, the A-staralgorithm (Hart et al., 1968), a modified method of Dijkstra’s algorithm,
uses heuristics by introducing a cost function for determining theshortest path to the destination (Fig. 1b). As compared to the original
Dijkstra’s algorithm, the A-star algorithm enables us to significantly
reduce computational cost to determine the optimal route by preferen­
tially searching for the adjacent links and nodes pointing to the desti­
nation. As illustrated above, there are a substantial number of studies
and algorithms for WR with some comparative assessments among themethods (e.g., Roh, 2013; Walther et al., 2016), but the definitive al­
gorithm suitable for WR has not yet been confirmed.
In this study, our first objective is to propose a new marine WR
system based upon the A-star algorithm for determining the shortest
path to the destination. A cost function is additionally introduced toaccount for oceanic and atmospheric conditions around the vessel ofinterest, combined with the heuristic function defined in the original Astar algorithm. In other words, this study proposes to modify the costfunction so that three options are considered to search for the optimal
paths: shortest travel distance, shortest travel time, and minimal fuelconsumption.
Furthermore, most of previous WR studies have focused on oceangoing vessels and ocean liners that cruise mostly on the open ocean.
Their routing algorithms are mainly for long voyages where a wide
range of route choices is possible. Nevertheless, it is difficult to apply

1. Develop a new WR system based on the A-star algorithm, whichcalculates an optimal route with consideration of oceanic and at­mospheric conditions, having capabilities of rerouting and avoidance
functions for more accurate and safer navigation;
2. Assess the efficiency of the newly developed WR system by applying
it to several scenarios including evacuation from a migrating
typhoon; and3. Configure a compact regional ocean circulation model executable on
typical laptop PCs aiming at a vessel-borne WR.
This paper is organized as follows. Section 2 describes the proposed
algorithms for the WR system. The configurations of the numericalmodels are given in Section 3. The experimental results with the new WR
algorithm are presented and discussed in Section 4. Finally, summariesof this study are presented in Section 5.
2. Weather routing algorithms
2.1. A-star and Dijkstra’s algorithms
In this study, a new marine ship routing system is proposed based on
the A-star algorithm (Hart et al., 1968, Fig. 1b). This algorithm finds theshortest path to the destination using inter-node traveling cost (hereafterχ ) and a heuristic function (hereafter h). The heuristic function h isdefined as the shortest path from each node to the final destination node
determined a priori to minimize the cost χ to search neighboring nodes
for candidates comprising the optimal path. Hence, the A-star algorithm
enables us to significantly reduce computational time needed to deter­mine the optimal route by using h as the “guidepost.” If h is set at zero,
the A-star algorithm simply reverts back to Dijkstra’s algorithm(Fig. 1a).
Fig. 2 compares the efficiency of the Dijkstra’s algorithm and the Astar algorithm in searching for the shortest path from the origin to thedestination by illustrating the estimated shortest paths and the numberof searches at every grid point. Even if no obstacles exist on the possible
path (Fig. 2a and b), Dijkstra’s algorithm with h ¼ 0 searched all the gridFig. 1. Schematic diagrams of searching algorithms of (a) original Dijkstra’s
method, (b) original A-star method, and (c) proposed new weighted A-starmethod. The start node denoted by m at the center seeks the optimal route
among the eight neighboring candidate nodes n1 – n8 during the first iteration.
Lengths of the arrows denote relative magnitudes of the cost χ . The star marks
are the final destinations located in the far field with heuristics hm at m. The Astar algorithm in (b) may preferentially search only n2 and n3 because of theshortest heuristics of hn2 or hn3. The ambient current approaches from the east
to alter the cost χ in (c), either to shorten (red arrows) or to elongate (bluearrows) the effective distance. (For interpretation of the references to color inthis figure legend, the reader is referred to the Web version of this article).
2

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Fig. 2. Schematic plots of Dijkstra’s algorithm
and the A-star algorithm examples. The black
curves are the shortest routes determined by
Dijkstra’s algorithm (a and c) and the A-star al­
gorithm (b and d). The green and magenta
squares are the origin and destination, respec­tively. The circles indicate that the grid points
are searched during iteration with colors indi­
cating number of the searched steps. The blacksquares are obstacles excluded from the graphsearches. (For interpretation of the references tocolor in this figure legend, the reader is referredto the Web version of this article).

points between the origin and the destination (Fig. 2a). On the other
hand, the A-star algorithm successfully finds the shortest path without
any unnecessary searches (Fig. 2b) that ensures significant improvement
in computational efficiency. If we place obstacles representing break­waters of port facilities near the origin and destination, the two algo­rithms detect the same optimal routes, whereas the A-star algorithm
requires fewer computational steps and searched grid points than
Dijkstra’s algorithm (Fig. 2c and d). Therefore, the A-star algorithm
significantly reduces search cost to find the shortest path without losing
precision as compared to Dijkstra’s algorithm.

Thus, in the shortest travel distance method, the optimal route issearched without considering oceanic currents, so the inter-node trav­
eling cost function from node m to node n (χ n|m) is defined merely asχ njm ¼ L;

where L is the arc length between nodes n and m as computed by Eq. (1).
In the second modification, we consider a weight on the cost χ byreflecting changes of the vessel’s relative velocity to the local back­ground currents as changes of an effective inter-node traveling distance(Fig. 1c); i.e., we assume that the traveling distance is either shortenedor elongated due to background current velocity. If a vessel is sailing on
still water at Vo but the ocean moves at Vi due to currents, then the
resultant absolute vessel speed is VoþVi. Accordingly, χ n|m is modified as:

2.2. Three options useful for marine WR
The traveling cost χ and heuristics h inherent in the original A-staralgorithm are designed with Euclidean distances. The optimal route isdetermined by minimizing the sum of the cost χ and the heuristics h,such that the optimal route follows the paths of shortest Euclidean dis­
tance. In this study, we incorporate the great-circular distance (herein­
after, simply referred to as an arc length) to χ and h, and the effect ofoceanic current and other factors to χ . The first modification has usually
been done in other marine WR algorithms. The arc length L between two
locations at (latitude φ1, longitude λ1) and (φ2, λ2) on the globe with adiameter of R may be approximated as
L ¼ R acos½sinφ1 sinφ2 þ cosφ1 cosφ2 cosðλ1

λ2 Þ�:(2)

χ njm ¼

Vo
L:Vo þ Vi(3)

If the ambient current direction coincides with the vessel’s traveling
direction, i.e., VoþVi > Vo, then the weighted χ n|m < L, so that the navi­
gator may favor sailing the vessel in this direction for the shorter traveltime.
In contrast, if the underlying current is in the same direction as thevessel, then the vessel can decelerate its relative speed to Vo Vi to keep
the absolute vessel speed unchanged to reduce fuel consumption. Ac­
cording to previous studies (e.g., Lu et al., 2015; Bialystocki and(1)
3

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Konovessis, 2016), fuel consumption of vessels is quadratically propor­
tional to the relative vessel speed. Hence, we introduce another modi­
fication for the minimal fuel consumption method, in which
�2
�
V
Vi
χ njm ¼ L o
:
(4)
Vo

Table 1Model configurations.

As a result, three options are introduced to search for the optimum
paths with the shortest travel distance: (1) without considering oceanic
currents, (2) minimizing travel time while maintaining a constant pro­peller speed, and (3) minimizing fuel consumption by preserving a
constant vessel speed affected by oceanic currents while altering pro­peller speed. Although the first method corresponds to the original Astar algorithm, the latter two methods rely on the weighted cost function
χ modified by oceanic and atmospheric conditions.
In addition, a rerouting function is also implemented, as well as an
avoidance function, for arbitrary areas where there are waves, wind,chemicals, and any conditions that are undesirable for navigation. The
former is achieved simply by re-computation of the optimal route
whenever it is demanded upon updates of oceanic, atmospheric, orvessel conditions. The latter is viable by removing the nodes that thevessel should avoid from WR computations.

ModelsNWPO

Seto-2km

Seto-600m

Computational
periodGrid cells

1/1/2013 - 12/31/
2014
800�480 (�32vertical layers)
2.0 km1/1/2013 - 2/28/
2013
224�192 (�20
layers)2.0 km9/1/2006 - 11/9/
2015
800�480 (�32
layers)
600 m

120 s

120 s

30 s

JMA-MSM (hourly)

JMA-MSM (hourly)

COAMPS bulk
formulaJCOPE2 (daily)

COADS (monthly
clim.)JCOPE2 (daily)

JMA-MSM (twohourly)COAMPS bulk
formulaParent ROMS-L1(daily)JCOPE2 (10-day
averaged)
SRTM30þJEGG500
Japan River
Association(monthly
climatology)
None

None

None

SRTM30þJEGG500
None

SRTM30þJEGG500
Japan River
Association(monthly clim.)

TPXO7.0

TPXO7.0

Horizontal grid
resolution
Baroclinic time
step
Surface windstress
Surface flux
Boundary/
Initialcondition
T-S nudging
Topography
Major riverdischarges
Tide

3. Regional ocean circulation models

weak nudging inverse time scale of 1/20 d 1 is applied towards the
10-day averaged JCOPE2 3D temperature and salinity fields to assure
long-term reproducibility of the Kuroshio path. This model was exten­
sively validated to show good agreement with satellite and in situ ob­
servations in Tada et al. (2018). The NWPO model is used for testinggraph-search problems in consideration of the Kuroshio path andrerouting and avoidance functions under a typhoon condition in Sec­
tions 4.1 and 4.2.

The ocean model we use in this study is the Regional Oceanic
Modeling System (ROMS; Shchepetkin and McWilliams, 2005, 2008), a
terrain-following, 3D primitive equation ocean circulation model withthe hydrostatic and Boussinesq approximations. The following sub­sections describe the model configurations used for the WR experiments
to validate and evaluate the proposed algorithms.
3.1. Northwestern Pacific Ocean model
The first configuration used in this study is the Northwestern Pacific
Ocean downscaling ROMS model (hereinafter NWPO; Fig. 3 andTable 1) developed in Tada et al. (2018). The daily averaged JCOPE2
reanalysis with a lateral grid spacing of ~10 km is used as the initial andopen boundary conditions. The atmospheric conditions are imposed byusing the Grid-Point-Value MesoScale Model (GPV-MSM), ahigh-resolution atmospheric reanalysis product published by the Japan
Meteorological Agency (JMA). The NWPO model domain extends hor­izontally 1,840 km � 1,440 km with a lateral grid spacing of 2 km toencompass a sufficiently wide area that includes the Kuroshio path andthe track of Typhoon 201418 (Phanfone). Temperature-salinity nudging(T-S nudging; Marchesiello et al., 2003; Uchiyama et al., 2018a) with a

3.2. Compact ocean model
We configure a compact ocean model for the SIS (hereinafter Seto2km; Fig. 3) that is sufficiently suited to execute on laptop computersor even on tablets and smartphones for a vessel-borne, quasi-standalone
WR system. The model is again based on ROMS and designed toencompass the SIS and the Kuroshio with a horizontal grid spacing of
2 km. A one-way offline nesting is employed from the assimilative
JCOPE2 oceanic reanalysis. Ten tidal constituents from the TPXO7.0
global tidal reanalysis (Egbert et al., 1994; Egbert and Erofeeva, 2002)
are additionally imposed on the open boundaries to account for sub­stantial tidal currents in the SIS. Other numerical conditions are listed inFig. 3. Model domains with bathymetry in color.
(a) Black box: JCOPE2, red box: NWPO model,
and magenta box: Seto-2km model (also shownin b), (b) black box: Seto-600m model. The black
dotted line in (b) identifies the transect A-B for
the cross-sectional plots in Fig. 5. The cyan cir­cles depict the locations of tide gauge stations
used for the harmonic analysis in Fig. 7. (Forinterpretation of the references to color in thisfigure legend, the reader is referred to the Web
version of this article).

4

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Table 1. By reducing horizontal and vertical resolutions, omitting river
inflows, and simplifying the surface flux estimation without losing therequired accuracy, the model can compute 3D oceanic currents anddensity field for one calendar day in approximately 4 min using a typicallaptop PC. In the following analyses, the Seto-2km ROMS model runs fora two-month period including the initial spin-up of one month.

each model domain are given as monthly climatological freshwater masspoint sources based on the Annual Record of Rainfall and Discharge
Database issued by the Japan River Association. The Seto-600m model
was also carefully validated with observed surface temperature andsalinity in Kosako et al. (2016) as well as high-frequency tidal signals in
the SIS in Uchiyama et al. (2018b). Other numerical configurations are
the same as those of the Seto-2km model, as summarized in Table 1.
It is anticipated that the Seto-2km model would be degraded bydesign from its higher-resolution counterpart due to the simplified
physics and numerical methods for the sake of portability. The point tobe confirmed here is whether the Seto-2km is still reasonably precise for
WR or not. Fig. 4 shows the comparison of sea-surface variables between
the Seto-2km and Seto-600m models averaged over one month inFebruary 2013. The surface temperature agreement is exhibited be­
tween the two models, whereas the agreement is moderate for salinityand velocity. The surface salinity of the Seto-2km model is over­estimated, especially in the inner-bay areas, because the compact model

3.3. Validation of the compact ocean model
The results of the Seto-2km model are compared with those ofanother higher-resolution SIS model developed in Kosako et al. (2016)(hereinafter Seto-600m; Fig. 3 and Table 1). The Seto-600m model is in a
double-nested configuration based on ROMS with grid spacings of 2 km
and 600 m, embedded in the JCOPE2 reanalysis. The surface fluxes are
obtained using the bulk formula that has been implemented in the
Coupled Ocean/Atmosphere Mesoscale Prediction System (COAMPS;
Hodur, 1997). The freshwater discharges from all the major rivers in

Fig. 4. Spatial maps of temperatures (a and b), salinities (c and d), and velocity magnitudes (e and f) at surface averaged over one month in February 2013. Leftpanels (a, c, and e) are the Seto-2km model results, while right panels (b, d, and f) are the Seto-600m model results.
5

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

does not account for major rivers that results in underestimation of
freshwater input to the SIS (Fig. 4c and d). The mean surface velocity
distributions of the compact model are almost identical to those of the
Seto-600m, despite the modest underestimation in the offshore areaaround the Kuroshio path, where the highest velocity emerges (Fig. 4e
and f). Such velocity intensification is well known to be associated withgrid-size refinement, because the frontal structures around the Kuroshio
path are better resolved to enhance cross-frontal lateral buoyancy andpressure gradient forces, leading to geostrophic axial acceleration of theKuroshio. The vertical cross-sectional plots along the A-B line (see Fig. 3
for the location) of the same variables as those in Fig. 4 show that the
Seto-2km model successfully reproduces the stratification of tempera­ture and salinity near the Kuroshio front (Fig. 5a and b). In contrast, thestreamwise velocity along the Kuroshio has an overall similarity, butwith intensified streamwise velocity in the Seto-600m due to theenhanced frontal gradient (Fig. 5e and f). As anticipated, although themodel grid resolution and precise forcing are important in ensuring thedetailed reproducibility, the Seto-2km model has a reasonable accuracy
comparable to the high-resolution model, except for some limited re­gions such as near river mouths and the Kuroshio path.
Fig. 6 compares the ocean currents of the Seto-2km model and thosefrom in situ observation data obtained by the Ministry of Land, Infra­
structure, Transport, and Tourism. This data contains current direction
and velocity observed by upward-looking ADCPs deployed at the bot­toms. The panels show comparisons in three layers, surface, middle, andbottom, at three points shown in Fig. 6a at 12-h intervals for one-month
period in February 2013. Root mean square differences (RMSDs) for
modeled currents relative to the in situ data are calculated as

RMSD ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffin1X
ðfyi Þ2 ;
n i¼1 i(5)where n, f, and y are the number of observations, modeled forecast
value, and observed value, respectively. At both P1 and P2 in Fig. 6a, theSeto-2km model reproduces currents in the bottom layer better thanthose in the surface and middle layers (Fig. 6 b-g). The RMSDs at P3
shows that the Seto-2km model estimated the current well at each layer(Fig. 6 h-j). Consequently, even though the Seto-2km model is designed
as a compact ocean model, its reproducibility of depth-dependent hor­izontal current velocity is reasonably well, while we expect significant
improvement by using finer resolution and more accurate and sufficient
initial condition and forcing datasets.
Fig. 7 shows the result of a tidal harmonic analysis of surface ele­vations from the Seto-2km model and in situ tidal observations at seventide gauge stations of JMA at Kobe, Wakayama, Komatsushima, Kochi,
Matsuyama, Uwajima, and Saiki, shown by the cyan markers in Fig. 3b
for a one-month period in February 2013. We show the result for four
selected principal tidal constituents (M2, S2, K1, and O1) as scatterplotsof amplitudes and phase epochs. It is clearly demonstrated that theamplitudes and phases are reasonably well reproduced by the Seto-2km
model with the correlation coefficients of 0.944 and 0.999, respectively.
This result assures us that the present Seto-2km model adequately re­produces high-frequency tidal currents prominent in the SIS. From theresults demonstrated above, we would conclude that the compact Seto2km model is exploitable for ship routing problems in the SIS.

Fig. 5. Cross-sectional plots of temperatures (a and b), salinities (c and d), and streamwise velocities normal to the cross-section (e and f) averaged over one month inFebruary 2013 along the transect A-B shown by the black dotted line in Fig. 3b. Left panels (a, c, and e) are the Seto-2km model results, while right panels (b, d, and f)
are the Seto-600m model results.
6

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Fig. 6. Spatial map of observation points and time series of current velocity vectors for P1 (b–d), P2 (e–g), and P3 (h–j). Panels (b, e, f), (c, f, i), and (d, g, j) show thedata obtained in surface, middle, bottom layer, respectively. The black and magenta vectors show the velocity from observation data and the Seto-2km model,
respectively. The observation data is available at 12-h intervals for one month in February 2013. The RMSDs for each point and layer are represented in the top
middle box.

Fig. 7. Scatter plots of four principal tidal constituents (M2, S2, K1, and O1) derived from the harmonic analysis of hourly surface elevations from Seto-2km model andobservations at seven tide gauge locations shown in Fig. 3b. (a) Amplitudes and (b) phase epochs. R in each plot indicates the correlation coefficient.

4. Experimentsfunction is also not considered in this experiment for brevity andsimplicity. The imposed speed of the vessel is 10 kn unless otherwise
stated.
Fig. 8 shows the three optimal routes and associated calculation steps
at each of the searched grid points. The black curve in Fig. 8a represents
the shortest travel distance route based on the original A-star algorithm.
Because this route is simply determined without any oceanic influences,
it is merely a straight line connecting the origin and the destination
while avoiding islands. We notice that the eastward zonal travel occurs
after reaching the coastal area for shorter arcs at higher latitude. By
contrast, the shortest travel time route (Fig. 8b) partially utilizes theKuroshio current flowing mostly parallel to the vessel’s track. The vessel
does not take a roundabout route on the Kuroshio off Shikoku Island, butrather takes a straight course before it starts using the Kuroshio from the

4.1. The three new WR methods: a case study for an ocean liner
In this subsection, the optimal routes are evaluated by using the three
new routing strategies. A test problem considered here is WR of an ocean
liner from Okinawa to Tokyo (Fig. 3b). A key point of this application is
how the vessel utilizes the Kuroshio drifting northeastward (Fig. 8).
Because the Kuroshio is recognized to be intense enough at up to 3 m/s(~6 kn) to affect a vessel’s speed, some previous studies have tried toincorporate the current effects in WR (Chen et al., 2015; Chang et al.,2013). Hourly-averaged oceanic surface currents are provided from the
NWPO model, whereas we used the values of surface currents at the
departure time that were frozen during the voyage. The rerouting
7K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Fig. 8. Optimal routes determined by the newly developed WR algorithms (a, b, and c) and numbers of the searched steps (d, e, and f). Black curves in (a and d), (band e), and (c and f) are the shortest distance route, shortest time route, and minimum fuel consumption route, respectively. Colors in upper plots indicate thesnapshots of the surface velocity magnitude, while those in lower plots are numbers of the searched steps. The green and magenta squares are the origin (Okinawa)
and destination (Tokyo), respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article).

Kii Peninsula. On the other hand, the minimal fuel consumption route(Fig. 8c) is quite different from the other two routes. The vessel uses theKuroshio for a long distance as much as possible to minimize fuelconsumption.
The lower panels of Fig. 8 show that the shortest travel time (Fig. 8e)and minimal fuel consumption methods (Fig. 8f) require much fewer
searches and thus are computationally more efficient to find the optimalroute than the original A-star method (Fig. 8d), even though they ac­
count for influences of the Kuroshio in the evaluation of the cost χ . The
two weighted A-star searches suggest detours to utilize the Kuroshio
following the vessel to reduce χ , whereas the minimal fuel route favors
much greater detours. Apparently, the degrees to which the weighted
methods recommend utilization of the following currents depend on theratio of the vessel’s speed to the ambient current speed. The Kuroshio
current is at 2–3 m/s (4–6 kn) near the ocean surface, which is signifi­
cant as compared to the vessel’s speed set at 10 kn. If the vessel were on
the Kuroshio path during its entire voyage, the expected fuel con­sumption could be minimum on average. However, such a route requires
large detours that eventually lead to extra fuel consumption due tolonger travel time. The weighted cost functions optimize the balance
between these competing influences on the routing purposes. On theother hand, whereas the weighted searches are computationally lesscostly than the A-star search, the shortest travel time search takes more
search steps to finalize the optimization than the minimum fuel search,
particularly in the east of the Kii Peninsula. This result illustrates asensitivity of the optimization to reduce the travel time by considering abalance between the faster vessel’s absolute speed due to following theKuroshio and the longer travel distance due to detours.
The results of the three optimal path findings are summarized inTable 2. Needless to say, the shortest travel distance route is the shortest.
In the case of the shortest travel time method, although the distancetraveled is longer than the shortest distance route by ~20 km, the al­
gorithm successfully shortens the travel time by approximately 4 h,
thanks to the Kuroshio effect. In the case of the minimum fuel con­sumption method, whereas the traveled distance and time are increased
greatly, the fuel consumption is reduced by ~8.4% compared with theshortest distance route. Consequently, it is confirmed that the developed
WR algorithms successfully determine the optimal vessel tracks with thethree options. It is worth noting that the present method with theweighted cost χ is quite flexible and readily expandable to cases under
Table 2
Traveled time, traveled distance, and fuel efficiency from the results of the three
optimal routes shown in Fig. 8. Note that fuel efficiencies are normalized by thatof the shortest distance route.

Shortest Distance Route
Shortest Time
Route
Minimum Fuel
Consumption Route8

Traveled Time(min)

Traveled
Distance (km)Fuel
Efficiency

5066.94812.7

1626.81647.1

1.000
0.969

5232.21732.30.916

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

other restrictions such as with retardation of the vessel’s speed due tosurface waves, winds, and even surrounding traffic.

dotted curve is the shortest travel time route without considering thetyphoon, while the blue dotted curves are the predicted shortest timeroutes sequentially updated with hourly rerouting. The black solid
curves are the actual vessel tracks. The optimal vessel track is pro­nouncedly modified with the updated weather information to success­fully avoid encountering the migrating typhoon. Interestingly, theresultant optimal route is on the right side of the typhoon where wind
and waves are usually larger than those on the other side in the northernhemisphere. Obviously, this optimal path is affected by the along-coast
typhoon track attributable to less severe weather conditions on the rightside of the typhoon for the vessel of interest. The vessel navigates to­wards Okinawa, not only keeping away from the typhoon but also
minimizing excess distance traveled to bypass the storm. These results
confirm that the present WR algorithm can determine desirable optimalroutes with the rerouting and avoidance functionalities.

4.2. Avoidance of a typhoon with the rerouting function
The rerouting and avoidance functions are tested for a case with atyphoon traveling near the vessel track sailing from Tokyo to Okinawa.
The same NWPO model used in the preceding application provides theoceanic conditions. An analysis is conducted with the shortest travel
time method for the period when Typhoon 201418 (Phanfone) emerges
in the domain. Many previous studies have tried to navigate ships inwindy and high-wave conditions by estimating the impact of the con­ditions on the vessel’s body and motion (Journ�ee and Meijers, 1980;
Shaoze et al., 2016). However, in this subsection, we account for suchsevere weather conditions as simply as possible to evade complexity
associated with interactions between the ship’s body and the ocean
merely for testing of the implemented avoidance function. The avoid­ance condition is set with two criteria, viz., wind speed greater than
17.2 m/s and wave height higher than 7.5 m. These weather conditions
are provided from the GPV-MSM atmospheric product and the
GPV-CWM spectral wave product, both published by the JMA. The
nodes that meet either of the two criteria are configured as “obstacles,”similar to the land nodes, where the vessel is unnavigable. The
time-varying ambient currents are also considered in this experiment.
The rerouting is also performed at 1 h intervals to synchronize with thehourly update of the weather and oceanic information.
Fig. 9 shows the snapshots of the migrating typhoon visualized bysignificant wave height (colors) and wind speed at 10 m above the ocean
surface (contours) with the vessel’s tracks at intervals of 8 h. The black

4.3. Ship routing under strong tidal currents
In this subsection, the Seto-2km compact ocean model is exploited todevelop and test a quasi-standalone WR system suitable for regional shipnavigation in estuaries, bays, and harbors. As argued in Section 3, thismodel is sufficiently capable of running on laptop computers, tablets,
and smartphones on the ship while maintaining reasonable accuracy.
Therefore, our motivation to employ such a compact model totally dif­fers from that in the preceding subsections, where the high-resolution
oceanic downscaling products originally designed for research pur­poses were used for precise evaluation of the proposed algorithms.
The SIS is characterized by complex topography that forms a typicalbasin-strait system (Fig. 3), where strong semi-diurnal tidal currents

Fig. 9. Snapshots of the migrating typhoon and optimal routes with the rerouting function at intervals of 8 h. Colors indicate significant wave heights from GPVCWM, and contours are surface wind speeds from GPV-MSM. The dotted black curves are the shortest travel time routes without considering the typhoon, while
the dotted blue curves are the predicted shortest time routes sequentially updated with hourly rerouting. The solid black curves are the actual vessel tracks. The green
and magenta squares are the origin (Tokyo) and destination (Okinawa), respectively. Green circle is the position of the vessel at the labeled time. (For interpretationof the references to color in this figure legend, the reader is referred to the Web version of this article).
9

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

develop near topographies (e.g., Uchiyama et al., 2018b; Zhang et al.,2019). We intend to find the optimal southwestward routes from Osaka
to Kagawa, which requires passing a narrow tidal strait where tidal
currents are quite energetic. Because tidal currents fluctuate largely at
the semi-diurnal period in the area, travel time and fuel efficiency areaffected by their direction and intensity that depend on tidal phase andthus on the departure time.
Fig. 10 shows the optimal routes from Osaka to Kagawa with theshortest travel time method for the ship operated at 10 kn leaving Osaka
every 6 h from 0:00 UTC on February 22, 2013. The color represents thehorizontal velocity magnitude averaged over the 6 h after each depar­
ture. The vectors correspond to the direction and speed of the timeaveraged tidal currents. Overall, the determined shortest travel timeroutes are similar to each other among all the cases. The present WR
system suggests passing through the northern strait with rather straightroutes from the origin to the strait and from the strait to the destination.
The most striking differences associated with the tidal phase areconfirmed by comparing two cases, in which the ship departs at 12:00
and 18:00 on February 22 in 2013 (Fig. 10c and d). As labeled in Fig. 10,
although the estimated routes, travel distances, and durations areapproximately the same, the difference of fuel efficiency is 14.5% for
these two cases. Fig. 11 shows the time series of the fuel efficiency ac­
cording to the departure times shown on the abscissa, normalized bythat with the least efficiency. The estimated fuel efficiency with different
departure times varies at a ~12-h period. We found that the overall fuelefficiency depends on the condition of the tidal currents at the northernstrait, rather than on the departure time. Apparently, this is caused by
the predominance of the semi-diurnal tidal variability in the SIS. For
instance, Fig. 10c clearly demonstrates that strong eastward tidal cur­rents occur at the strait, where the narrow channel topography in­tensifies tidal currents opposing the westward traveling vessel, resulting
in deceleration of the vessel. In turn, if the vessel leaves the origin at
18:00 (Fig. 10d), it approaches the strait at ~20:00, when the tidal
currents are headed westward to accelerate the vessel. From theseFig. 10. Shortest time routes from Osaka (the origin, green square marks) to Kagawa (the destination, magenta square marks) with different departure times
occurring every 6 h from 0:00 UTC on February 22, 2013. Colors indicate the surface velocity magnitudes averaged over 6 h after each departure. The black vectors
are superposed to show the subsampled time-averaged surface current velocity. (For interpretation of the references to color in this figure legend, the reader isreferred to the Web version of this article).
10

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706

Fig. 11. Fuel efficiencies for labeled departure times in the abscissa, starting at 0:00 UTC on February 22, 2013. The fuel efficiency is relative to that at 17:00 UTC(red circle). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article).

for arbitrary oceanic and weather conditions around the vessel. Possible
modifications are to be made by adjusting the weight defined as afunction of the vessel’s speed. For instance, given detailed wave or windconditions around the vessel, the traveling speed may be alteredaccordingly, which in turn influences the effective inter-node distanceequivalent to the weighted cost. The challenge towards an integrated
vessel-borne WR system with a compact ocean model is also a uniqueundertaking of this study. Such a system allows us tremendous flexibility
in utility of the marine WR, even though data retrieval is impossible orlargely limited if the vessel is out of radio range under certain circum­
stances, including cases of severe natural disasters.

results, the direction of the semi-diurnal tidal currents in such a small
estuary has a pronounced impact on vessel speed, leading to substantial
alteration of fuel efficiency.
5. Conclusions
This study aimed to develop a weather routing (WR) system forvessel navigation based on the A-star algorithm. Building upon theheuristics introduced to reduce computational cost as compared to thatof the widely used Dijkstra’s algorithm, the cost function was addi­
tionally modified to account for oceanic and atmospheric conditions
around the vessel of interest. Three options were implemented to search
for the optimal paths with the shortest travel time, minimal fuel con­sumption, and shortest travel distance. In addition, we introduced arerouting function that recalculates and updates the optimal routerecursively whenever the route should be reconsidered with updated
oceanic and weather information. An avoidance function to exclude
specific nodes from the search was also incorporated into the WR system
by setting thresholds according to oceanic and weather conditions. The
ocean liner examples along the northeastward drifting Kuroshio current
off Japan with and without a migrating typhoon were chosen to test thedeveloped WR system and to confirm the expected functionalities.
Furthermore, a compact regional ocean circulation model, which is
sufficiently accurate and executable on typical laptop PCs, was config­
ured for a vessel-borne WR system. This compact model was successfully
applied to evaluate the optimal vessel paths in the SIS, where highfrequency tidal currents modified by complex topographies are signifi­
cant enough to alter the vessel’s absolute speed, to which fuel efficiency
is eventually attributed to a remarkable extent. In this study, our pur­pose was accomplished by using simple datasets to design the compact
ocean model. For further applying this method to the practical imple­
mentation of WR systems in the Japan’s neighboring seas in the future,
the compact model will need to be operational and also coupled withatmospheric models. In such a case, the compact model would become
more complicated and time-consuming, because forecasts of winds,
waves, and currents for several days up front are necessary. In otherwords, more accurate and sophisticated WR system will be realized withhigher-performance PC and more frequent data updates.
The developed WR system is quite general, computationally very low
cost, and universally applicable to any vessel navigation scenarios, notonly in the open ocean but also in estuaries and harbors with complex
topographies. Although we mostly examined the effects of the ambient
surface currents, the proposed weighted-cost function based on the Astar algorithm has substantial flexibility and expandability to accountAcknowledgments
This study was financially supported by the Japan Society for the
Promotion of Science (JSPS) Grants-in-Aid for Scientific Research
15KK0207 and 18H03798 at Kobe University. We also appreciate Yota
Suzue (CTI Engineering Co., Ltd.), Hiroaki Tada (Nippon Koei Co. Ltd),and Yasumasa Miyazawa (Japan Agency for Marine-Earth Science andTechnology) for their help in the development of the ROMS models
embedded in the assimilative Japan Coastal Ocean Predictability Ex­
periments (JCOPE2) oceanic reanalysis that was used in the present
study.
References
Andersson, A., 2015. Multi-objective Optimisation of Ship Routes. Master’s thesis.
Chalmers University of Technology, G€
oteborg, Sweden.
Bialystocki, N., Konovessis, D., 2016. On the estimation of ship’s fuel consumption andspeed curve: a statistical approach. J. Ocean Eng. Sci. 1 (2), 157–166. https://doi.
org/10.1016/j.joes.2016.02.001.
Bouman, E.A., Lindstad, E., Rialland, A.I., Strømman, A.H., 2017. State-of-the-art
technologies, measures, and potential for reducing GHG emissions from shipping – areview. Transp. Res. D Transp. Environ. 52 (A), 408–421. https://doi.org/10.1016/j.
trd.2017.03.022.
Buhaug, Ø., Corbett, J.J., Endresen, Ø., Eyring, V., Faber, J., Hanayama, S., Lee, D.S.,
Lee, D., Lindstad, H., Markowska, A.Z., Mjelde, A., Nelissen, D., Nilsen, J.,Pålsson, C., Winebrake, J.J., Wu, W.Q., Yoshida, K., 2009. Second IMO GHG Study
2009. International Maritime Organization, London, UK.
Calvert, S., Deakins, E., Motte, R., 1991. A dynamic system for fuel optimisation transocean. J. Navig. 44 (2), 233–265. https://doi.org/10.1017/S0373463300009978.
Chen, C., Shiotani, S., Sasa, K., 2013. Numerical ship navigation based on weather andocean simulation. Ocean Eng. 69, 44–53. https://doi.org/10.1016/j.
oceaneng.2013.05.019.
Chen, C., Shiotani, S., Sasa, K., 2015. Effect of ocean currents on ship navigation in theeast China sea. Ocean Eng. 104, 283–293. https://doi.org/10.1016/j.
oceaneng.2015.04.062.
Chang, Y.C., Tseng, R.S., Chen, G.Y., Chu, P.C., Shen, Y.T., 2013. Ship routing utilizing
strong ocean currents. J. Navig. 66, 825–835. https://doi.org/10.1017/
S0373463313000441.

11

K. Kurosawa et al.

Ocean Engineering 195 (2020) 106706
Roh, M.I., 2013. Determination of an economical shipping route considering the effects
of sea state for lower fuel consumption. Int. J. Nav. Archit. Ocean Eng. 5 (2),
246–262. https://doi.org/10.3744/JNAOE.2013.5.2.246.
Sen, D., Padhy, C., 2010. Development of a ship weather-routing algorithm for specificapplication in North Indian Ocean Region. Int. Conf. Mar. Technol. (MARTEC) 2010,21–27.
Shaoze, L., Ning, M., Hirakawa, Y., 2016. Evaluation of resistance increase and speed lossof a ship in wind and waves. J. Ocean Eng. Sci. 1 (3), 212–218. https://doi.org/
10.1016/j.joes.2016.04.001.
Shchepetkin, A.F., McWilliams, J.C., 2005. The regional ocean modeling system (ROMS):
a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean
Model. 9 (4), 347–404. https://doi.org/10.1016/j.ocemod.2004.08.002.
Shchepetkin, A.F., McWilliams, J.C., 2008. Computational kernel algorithms for finescale, multiprocess, longtime oceanic simulations. Handb. Numer. Anal. 14,
119–181. https://doi.org/10.1016/S1570-8659(08)01202-0.
Smierzchalski, R., Michalewicz, Z., 2000. Modeling of a ship trajectory in collisionsituations at sea by evolutionary algorithm. IEEE Trans. Evol. Comput. 4 (3),
227–241. https://doi.org/10.1109/4235.873234.
Szlapczynska, J., 2007. Multiobjective approach to weather routing. Trans. Nav. Int. J.
Mar. Navig. Saf. Sea Transp. 1 (3), 273–278.
Tada, H., Uchiyama, Y., Masunaga, E., 2018. Impacts of two super typhoons on the
Kuroshio and marginal seas on the Pacific coast of Japan. Deep. Res. Part I 132,
80–93. https://doi.org/10.1016/j.dsr.2017.12.007.
Trodden, D.G., Murphy, A.J., Pazouki, K., Sargeant, J., 2015. Fuel usage data analysis for
efficient shipping operations. Ocean Eng. 110 (B), 75–84. https://doi.org/10.1016/j.
oceaneng.2015.09.028.
Tsou, M.-C., Cheng, H.-C., 2013. An ant colony algorithm for efficient ship routing. Pol.
Marit. Res. 20 (3), 28–38. https://doi.org/10.2478/pomr-2013-0032.
Uchiyama, Y., Kanki, R., Takano, A., Yamazaki, H., Miyazawa, Y., 2018. Mesoscale reproducibility in regional ocean modeling with a 3-D stratification estimate based on
Aviso-Argo data. Atmos.-Ocean 56 (4), 212–229. https://doi.org/10.1080/
07055900.2017.1399858.
Uchiyama, Y., Zhang, X., Suzue, Y., Kosako, T., Miyazawa, Y., Nakayama, A., 2018.
Residual effects of treated effluent diversion on a seaweed farm in a tidal strait using
a multi-nested high-resolution 3-D circulation-dispersal model. Mar. Pollut. Bull.
130, 40–54. https://doi.org/10.1016/j.marpolbul.2018.03.007.
Walther, L., Rizvanolli, A., Wendebourg, M., Jahn, C., 2016. Modeling and optimization
algorithms in ship weather routing. Int. J. e-Navigation Marit. Econ. 4, 31–45.
https://doi.org/10.1016/j.enavi.2016.06.004.
Wei, S., Zhou, P., 2012. Development of a 3D dynamic programming method for weather
routing. Trans. Nav. Int. J. Mar. Navig. Saf. Sea Transp. 6 (1), 79–85.
Zhang, X., Uchiyama, Y., Nakayama, A., 2019. On relaxation of the influences of treatedsewage effluent on an adjacent seaweed farm in a tidal strait. Mar. Pollut. Bull. 144,
265–274. https://doi.org/10.1016/j.marpolbul.2019.04.050.

De Wit, C., 1990. Proposal for low cost ocean weather routing. J. Navig. 43 (3), 428–439.
https://doi.org/10.1017/S0373463300014053.
Dijkstra, E.W., 1959. A note on two problems in connexion with graphs. Numer. Math. 1(1), 269–271. https://doi.org/10.1007/BF01386390.
Efentakis, A., Pfoser, D., Voisard, A., 2011. Efficient data management in support ofshortest-path computation. In: Proc. The 4th ACM SIGSPATIAL International
Workshop on Computational Transportation Science. ACM Press, pp. 28–33. https://doi.org/10.1145/2068984.2068990.
Egbert, G.D., Erofeeva, S.Y., 2002. Efficient inverse modeling of barotropic ocean tides.
J. Atmos. Ocean. Technol. 19 (2), 183–204.
Egbert, G.D., Bennett, A.F., Foreman, M.G.G., 1994. TOPEX/POSEIDON tides estimatedusing a global inverse model. J. Geophys. Res. 99 (C12), 24821–24852. https://doi.
org/10.1029/94JC01894.
Hagiwara, H., 1989. Weather Routing of Sail-Assisted Motor Vessels. Ph. D. thesis. Delft
University of Technology, The Netherlands.
Hodur, R.M., 1997. The naval research laboratory’s coupled ocean/atmospheremesoscale prediction system (COAMPS). Mon. Weather Rev. 125, 1414–1430.
Hart, P.E., Nilsson, N.J., Raphael, B., 1968. A formal basis for the heuristic determinationof minimum cost paths. IEEE Trans. Syst. Sci. Cybern. 4 (2), 100–107. https://doi.
org/10.1109/TSSC.1968.300136.
Journ�ee, J.M.J., Meijers, J.H.C., 1980. Ship routing for optimum performance: Part I
Prediction of speed and power of a ship in a seaway. Trans. IME 12–21.
Jung, J.S., Rhyu, K.S., 1999. A study on the optimal navigation route decision using A*
algorithm. J. Korean Inst. Off. Autom. 4 (1), 38–46.
Kim, K., Lee, K.M., 2018. Dynamic programming-based vessel speed adjustment forenergy saving and emission reduction. Energies 11 (5), 1273. https://doi.org/
10.3390/en11051273.
Kosako, T., Uchiyama, Y., Mitarai, S., 2016. A multi-year analysis on circulation andassociated larval transport in the Seto Inland Sea, Japan. J. Jpn Soc. Civil Eng. Ser.
B2 72 (2), 1273–1278 (in Japanese with English abstract).
Kosmas, O.T., Vlachos, D.S., 2012. Simulated annealing for optimal ship routing.
Comput. Oper. Res. 39 (3), 576–581. https://doi.org/10.1016/j.cor.2011.05.010.
Lu, R., Turan, O., Boulougouris, E., Banks, C., Incecik, A., 2015. A semi-empirical ship
operational performance prediction model for voyage optimization towards energy
efficient shipping. Ocean Eng. 110 (B), 18–28. https://doi.org/10.1016/j.
oceaneng.2015.07.042.
Marchesiello, P., McWilliams, J.C., Shchepetkin, A.F., 2003. Equilibrium structure anddynamics of the California current system. J. Phys. Oceanogr. 33, 753–783.
Miyazawa, Y., Zhang, R., Guo, X., Tamura, H., Ambe, D., Lee, J., Okuno, A.,
Yoshinari, H., Setou, T., Komatsu, K., 2009. Water mass variability in the Western
North Pacific detected in 15-year eddy resolving ocean reanalysis. J. Oceanogr. 65(6), 737–756. https://doi.org/10.1007/s10872-009-0063-3.
Papanikolaou, A., Zaraphonitis, G., Bitner-Gregersen, E., Shigunov, V., El Moctar, O.,
Guedes Soares, C., Reddy, D.N., Sprenger, F., 2016. Energy efficient safe SHip
OPERAtion (SHOPERA). Transp. Res. Procedia. 14, 820–829. https://doi.org/
10.1016/j.trpro.2016.05.030.

12