# Long-Range Route-planning for Autonomous Vehicles in the Polar Oceans

Maria Fox1, Michael Meredith1, J. Alexander Brearley1, Dan Jones1 and Derek
Long2

(1British Antarctic Survey
2Schlumberger Cambridge Research
)

###### Abstract

There is an increasing demand for piloted autonomous underwater vehicles
(AUVs) to operate in polar ice conditions. At present, AUVs are deployed from
ships and directly human-piloted in these regions, entailing a high carbon
cost and limiting the scope of operations. A key requirement for long-term
autonomous missions is a long-range route planning capability that is aware of
the changing ice conditions. In this paper we address the problem of
automating long-range route-planning for AUVs operating in the Southern Ocean.
We present the route-planning method and results showing that efficient, ice-
avoiding, long-distance traverses can be planned.

This paper is under review111This work has been submitted to the Journal of
Atmospheric and Oceanic Technology. Copyright in this work may be transferred
without further notice..

##  1 Introduction

In the last decade, the rising demand for oceanographic data, to inform
diverse contemporary issues such as climate change, ecosystems services and
marine conservation, has spurred the development of autonomous platforms for
sustained acquisition of data. To extend the capability of existing networks
such as Argo [10], there is considerable interest in increasing the
coordinated use of piloted autonomous underwater vehicles (AUVs) [12] in the
hostile conditions of the polar oceans. Such vehicles, including low-power
underwater gliders, offer distinct advantages to ocean observing networks.
They can be piloted remotely and hence can conduct targeted campaigns in
specific areas, with tasking of the vehicles changing during missions in
response to evolving environmental conditions or user requirements.

At present, such autonomous underwater vehicles are commonly used individually
or in small clusters, with direct human piloting of individual vehicles from
land. There is strong interest in progressing capability to the point where
these vehicles can be deployed and controlled fully autonomously, in complex
environmental conditions and over long periods of time. For this capability to
be achieved in the polar oceans, the path-planning methods used must take
account of the changing ice conditions. In this paper we address the problem
of automating long-range route-planning for AUVs operating in the Southern
Ocean.

##  2 Background: Automated Long-distance Route-planning

We address the problem of finding a navigable path between two locations,
whilst respecting a collection of spatial and temporal environment
constraints, and optimising an objective function.

There is a large body of work in automated path planning for robotic vehicles,
including underwater vehicles. The most common methods use grid-based
discretisations or continuous numerical approaches. Grid-based methods
discretise the environment into a graph of locations and then use optimisation
methods to find paths from the start node to the destination. Approaches to
the problem of finding shortest paths in grids use heuristic search [2] or
greedy methods such as Dijkstra’s algorithm [3, 11]. Using a graph-based
representation, in which the environment is gridded at some chosen resolution,
edges can be weighted to model environmental impacts on travel-time, such as
the presence of strong surface currents. A shortcoming of grid-based
approaches is that all paths follow the edges of the pre-defined grid, between
vertices that are the centres and corners of the cells. This results in so-
called Manhattan paths, consisting of straight horizontal, diagonal and
vertical lines between cell centres, which are often far from optimal
depending on the grid resolution. Further increasing grid resolution to better
approximate optimal solutions is often prohibitively computationally
expensive. A further difficulty arises when environmental constraints are
time-dependent. These cannot easily be combined into edge weights, because
they change over time whilst the edge weights are constant. The 2D grid used
by Dijkstra’s algorithm must be extended to a 3rd dimension to model temporal
variability, which introduces extra complexity.

An alternative, cell-free, numerical modelling approach is to use the
principle of wave propagation to identify the points in an inhomogeneous space
that can be reached in equal time (or energy expenditure) from the location of
the vehicle. The approach avoids discretisation of the space at the outset,
but faces similar problems of efficiency and optimality of solutions that
depend on the number of points identified on each wavefront. Furthermore the
need for a temporal dimension must still be addressed when conditions are
time-dependent.

Both the grid-based and wave propagation methods have been used in long
distance route planning for both AUVs and ships [1, 8, 9, 6, 13]. Long-range
AUV operations in non-polar conditions are well-established: the Tethys Long-
Range AUV has successfully completed many science surveys consisting of long
traverses in the Pacific Ocean, of up to 3000km [5], using a path-planning
method based on Dijkstra’s algorithm [8] with the goal of optimising energy
efficiency.

Path-planning methods for underwater vehicles have thus far not considered
long-distance route planning in environments containing sea ice or icebergs.
Sea ice extent in the polar oceans varies month-by-month and year-by-year,
presenting many dangers and challenges to autonomous platforms. The historical
pattern of ice extent growth and retreat for the Antarctic is shown in Figure
1. Path-planning in these regions must take these seasonal patterns into
account as well as the behaviours of the ice in areas characterised by
perennial sea ice, such as the Weddell Sea.

Figure 1: The seasonal pattern of sea ice extent in the Antarctic. The blue
line is the mean sea ice extent over the years 1978-2021. Two standard
deviations are shown. The black line is 2013, a year of maximal ice extent.
The green line is 2010. The red line is 2016, a minimal year. The orange line
is 2021.

##  3 Long-distance Route-planning in Currents and Ice

Our objective is to automate long-distance route-planning for AUVs in the
polar oceans. AUV pilots typically focus on keeping well clear of the ice
edge, unless the AUV in question has specific ice-coping capabilities, which
most commercially available products do not. Our autonomously planned routes
will likewise avoid dangerous ice conditions whilst exploiting opportunities
to travel as close to the ice edge as is possible without compromising the
safety of the vehicle.

In the approach described in this paper we grid the environment and use
Dijkstra’s algorithm to find routes that minimise travel time. We address the
limitations of Dijkstra’s algorithm, noted above, in three ways. First, we
support a non-uniform gridding so that cells can be refined at the ice edge.
Second, we mitigate the Manhattan path problem by using Newton’s method to
find the optimal crossing point on the edge between two adjacent cells being
passed through on a route. This effectively allows an infinite set of vertices
along every edge, enabling substantial smoothing of routes without significant
computational expense (Newton’s method quickly converges to the optimal
crossing point on the edge). Finally, we add a temporal dimension by
discretising time into monthly steps and recomputing the gridding and
optimisation steps in each month. This combination leads to a computationally
inexpensive 3D (two spatial dimensions and one temporal dimension) shortest
path method that gives good estimates of the costs implied by AUV journeys
between distant waypoints. These can be readily recomputed during operation,
in response to unexpected changes in conditions, giving an adaptive route-
planning capability.

In this paper we focus on refining the grid along the ice edge, but the
refinement principle can also be used to refine the grid in other contexts,
for example: in the presence of land, ocean fronts, eddies, sharp bathymetric
gradients and so on. A principle of homogeneity is defined, which determines
the subdivision process. Cell sub-division is performed whenever a cell (of
any size) is determined to be inhomogeneous with respect to a specific
characteristic of interest. Routes can then be planned using Newton’s method
to find paths between adjacent cells and Dijkstra’s algorithm to find optimal
routes between waypoints through the non-uniform grid. We refer to this
combination of Dijkstra’s and Newton’s methods as the 3D dynamic shortest path
(3D-DSP) method.

The cost function may be any metric objective, such as travel time, distance,
power demands, fuel consumption or others, or indeed some combination of
metrics. In the work presented here, we optimise travel time, as this is an
important consideration for route-planning in dynamic conditions, and can also
be interpreted as a simple proxy for energy consumption. The 3D-DSP method can
be used in the spatially and temporally varying conditions of the polar
oceans, to calculate the fastest paths and travel-times, for autonomous
vehicles such as gliders, powered AUVs and other controllable vehicles that
may form part of a fleet.

##  4 Environmental Setting, Data and Methods

In this work, we use the Southern Ocean as an example domain. The region of
interest extends from 130∘W to 30∘E, and 40∘S to 80∘S (i.e. the Atlantic
sector of the Southern Ocean including the Weddell Sea). This region, an area
of more than 2.8 million km2, was chosen since it contains all of the natural
hazards present in any area of the global ocean, plus additional factors such
as sea ice. The model of currents and ice extent is built using
observationally constrained and physically consistent current data and sea ice
concentration data from the Biogeochemical Southern Ocean State Estimate
(B-SOSE) reanalysis product [7, 14].

Figure 2: Map showing the waypoints (blue numbered dots) of scientific
interest listed and named in Appendix B. The unlabelled red points are
scattered additional waypoints.

###  4.1 Representation of the ocean environment

The model of the environment is described in three stages. First, we describe
an initial procedure for coarsely gridding the region of interest shown in
Figure 2. The second stage is to introduce a temporal abstraction on top of
the spatial abstraction of the initial coarse grid. This allows the
representation of features that change over time, such as surface ice cover.
The third stage is to identify areas within the spatio-temporal abstraction
where grid refinement is needed in order to capture a description of the local
sea ice conditions that is adequate for safe route-planning.

####  4.1.1 Initial coarse grid

The B-SOSE data set used for this example contains current vectors each at a
specified coordinate in the Southern Ocean. We start by constructing an
abstraction of the field into a 2-dimensional grid of cells, each one
associated with an array containing the vectors belonging to that cell. To
characterise the cell, we compute its average current vector. Each cell
contains a centre point which is used for route-planning. Figure 3 shows an
arrangement of cells with the averaged current vector in each cell.

Figure 3: Part of the coarse grid abstraction of the vector field. Each cell
contains one vector which is an average of the currents in that cell.

Cell dimensions of 5∘ longitude by 2.5∘ latitude are use (this is
configurable). Given that the input file contains vectors at intervals of one
sixth of a degree, the maximum number of vectors that a 5 $\times$ 2.5 degree
cell can contain is 450. If it contains many fewer than this, then the cell is
likely to lie at least partially on land and may therefore not be navigable.
In the coarse grid a cell is classified as land if it contains fewer than 25%
of the maximum number of vectors.

####  4.1.2 Spatio-temporal abstraction

For each year in our data set, we extend the spatial abstraction described in
4.1.1 to take into account the presence and extent of sea ice in the same
region. The ice cover each year has a strong seasonal pattern, as is clear in
Figure 1, although this pattern has been changing over the last decade. For
route-planning, a monthly abstraction has been chosen. This is detailed enough
to enable the route planner to take into account the seasonal pattern without
becoming overwhelmed by detail on other timescales.

Each cell in the initial coarse grid, as well as recording the average current
vector, now records a list of 12 monthly abstractions. A monthly abstraction
contains the mean and variance of the sea ice area at that location over the
corresponding month. This extends the 2-dimensional representation into a
3-dimensional one in which there are twelve layers, each a grid representing
the region of interest in a given month.

A simple conservative threshold ice cover is used for determining that a cell
is ice-locked and therefore inaccessible. Whether or not a cell is ice-locked
changes from month to month. This means that some cells are navigable during
certain months and not others, depending on the year being modelled.

####  4.1.3 Refining the grid

A cell of size 5 $\times$ 2.5 degrees covers an area of more than 154,000km2
at the Equator. Cells closer to the pole are much smaller than this, for
example a cell at 60∘S covers about 77,000 km2. This is still extremely large
and adequate only for describing huge, deep, invariant regions of open sea
where the variance in current is small and there is no ice present. A non-
uniform grid is proposed so that, in complex regions such as the Weddell Sea,
cells of the grid are recursively subdivided.

We assume that refinement of the grid is based on the presence of data that
characterises some property or properties of each of the cells. Where that
data can be adequately characterised by a single model within a cell, we do
not need to refine it further, but in the cases where the data cannot be
appropriately characterised, the cell is decomposed. The spatial decomposition
is based on the use of a quad-tree [4] (widely used in representation of
2-dimensional data, such as images), refining a cell into quarters. The data
in the cell is then reallocated amongst the subdivisions and the process of
subdivision applied recursively to each of the quarters. We refer to the
requirement that must be satisfied in order to avoid further need for
subdivision as a homogeneity condition – homogeneous data needs no further
division, while inhomogeneous data requires that the containing cell be split.
This process can be continued arbitrarily far, but, in practice, once the data
has been split a number of times, the data content in each cell can become too
sparse to make further subdivision practical. Therefore, the definition of
homogeneity of the data typically includes a threshold on the quantity of
data.

Definition 1 shows how we define a homogeneity condition, $H$, for our
context.

###### Definition 1

A cell is considered homogeneous if, for every data set $D$ recorded in the
cell, a statistical measure, $g(D)$, is below the homogeneity threshold.

The splitting process is implemented using Algorithm 1. The input to the
algorithm is:

  * •

a cell, $C$, with bounds $a$ and $b$ (latitudes) in the $x$-dimension, and $c$
and $d$ (longitudes) in the $y$ dimension;

  * •

a data set of values $D$, with each datum associated with a point in $C$ by
the function $f$;

  * •

a boolean function, $H$, called the homogeneity condition;

$H$ is a boolean function that decides whether or not a cell should be split.
The statistical measure required to compute $H$, and the associated threshold,
will vary depending on the feature being considered. A detailed example is
given in Section 4.1.4, in which $H$ returns true (the cell is homogeneous and
should not be split) if the data set of ice concentrations in the cell falls
outside a specified range. If the homogeneity condition is not satisfied, then
the cell is split.

Algorithm: Split(C,D,f,H)

Data: A 2-dimensional cell, $C=\\{(x,y)\,|\,a\leq x<b,c\leq y<d\\}$, a
dataset, $D$, a function, $f:D\rightarrow C$, associating data in $D$ with
points in $C$ and a homogeneity condition, $H:\mathcal{P}(D)\rightarrow
Boolean$ on subsets of data in $D$

Result: A partition of $C$ into cells $\mathcal{C}$

begin

$\mathcal{C}\longleftarrow\emptyset$;

if  _H(D)_ then

add $C$ to $\mathcal{C}$;

return _$\mathcal{C}$_ ;

else

$C_{1}\longleftarrow\\{(x,y)\,|\,(x,y)\in C,x<(a+b)/2,y<(c+d)/2\\}$;

$C_{2}\longleftarrow\\{(x,y)\,|\,(x,y)\in C,x<(a+b)/2,y\geq(c+d)/2\\}$;

$C_{3}\longleftarrow\\{(x,y)\,|\,(x,y)\in C,x\geq(a+b)/2,y<(c+d)/2\\}$;

$C_{4}\longleftarrow\\{(x,y)\,|\,(x,y)\in C,x\geq(a+b)/2,y\geq(c+d)/2\\}$;

for  _i = 1 to 4_ do

$D_{i}\longleftarrow\\{d\,|\,d\in D,f(d)\in C_{i}\\}$;

$\mathcal{C}\longleftarrow\mathcal{C}\cup Split(C_{i},D_{i},f,H)$;

return _$\mathcal{C}$_ ;

Algorithm 1 Subdivide Cell, $C$ into partitions according to homogeneity
condition, $H$, applied to content data $D$.

Note that Algorithm 1 is recursive, allowing each subcell to be independently
further subdivided. The condition, $H$, checks both that there is sufficient
data in the subcell to support continued subdivision and also that the data
satisfies the required statistical property.

The construction of a non-uniform grid based on the homogeneity principle
allows us to combat the deficiencies of the coarse grid, while retaining the
benefits of a sparse structure where possible.

####  4.1.4 Inhomogeneity in the local ice concentration

We define a range, from a lower bound $lb$ to an upper bound $ub$, and a
threshold, $t$. Then, a cell month is considered inhomogeneous if between $lb$
and $ub$ of the ice measurements in that cell month are at $t$% or higher. If
the proportion of ice in the cell month above the $t$% concentration is below
$lb$%, we consider the cell month to be homogeneous open water: such a cell
can be navigated through in that month. At the other end of the range, if the
proportion is greater than $ub$%, then the cell month is considered
homogeneous ice. If the proportion is between these bounds, then the cell is
inhomogeneous and must be split so that the homogeneous sub-cells can be
found.

Given a cell month, $C$, and a data set of sea ice concentrations, $D$,
associated with $C$, we define a statistical measure, $g(D)$, as follows:

$g(D)=(lb-x)(x-ub),\,\,\,{\mbox{w}here}\,\,x=\frac{|\\{d\in
D\,|\,ice(d)>t\\}|}{|D|}$ (1)

The function $ice(d)$ returns the sea ice concentration (as a percentage)
associated with the datum $d$. The homogeneity condition, $H(D)$, which
decides whether the cell should be split, is $g(D)>0$, which is true whenever
the proportion of ice above the threshold lies inside the specified range. The
sub-splitting method quickly focuses on the parts of the cell where the high
percentages of ice reside in the cell and month being considered, enabling the
vehicle to traverse the ice edge more nimbly than is possible with the uniform
grid abstraction.

The degree of refinement of the grid impacts on waypoint accessibility and
degree of exposure to risk, therefore we consider two different homogeneity
conditions: Weak homogeneity: $lb=15$, $ub=70$, and Strong homogeneity:
$lb=5$, $ub=90$. A $4\%$ threshold is used in both cases, in order to focus on
the impact of the size of the range.

The intuition is that a homogeneity condition is weak if $g(D)>0$ is hard to
satisfy (the range is small), as this results in a lower level of
discrimination between different ice conditions. It is strong if, by contrast,
$g(D)>0$ is easy to satisfy (the range is large), as this encourages sub-
division and hence differentiation. When $lb=ub$, the homogeneity condition is
always true and the grid is never refined.

Since subdivision divides a cell into quarters, only 8 subdivisions are
required to reduce a 5 $\times$ 2.5 degree cell, of area 154,000 km2 at the
Equator, to an array of cells the smallest of which is 2 km2 in area. In
practice, the number of subdivisions possible is limited by the quantity of
data available in the resulting cells. Given that the B-SOSE data set supplies
measurements at one sixth of a degree, after 3 subdivisions a subcell contains
only 7 daily data points and, hence, 210 points for any given month. If
further splitting is performed, the ability to perform meaningful statistical
analyses is lost. We therefore restrict the depth of subdivision to 3 for the
experiments reported in this paper.

The degree of refinement of the grid varies between layers, reflecting the
fact that the sea ice area changes from month to month. Additionally, the sea
ice changes between years, so the 3-dimensional grids representing different
years may have very different levels of refinement in each of the layers.

###  4.2 The 3D-DSP Method

The grid abstraction represents the region as a graph,

$G=\langle V,E\rangle$ (2)

where $V$ is the set of vertices (cell centre points) and $E$ is the set of
edges (lines connecting adjacent pairs in $V$). As part of the problem
specification we are given a special set of coordinates called the waypoints
set (a complete list is provided in Appendix B). These are identified
locations with the special significance that the routes to be planned will
start and end at them. Every point in the waypoints set exists inside the cell
whose central point is closest to the waypoint coordinate.

The full 3D-DSP method involves two stages. First, the shortest paths are
computed from each source waypoint to all destination waypoints within each
month (the planning month), and hence each layer, of the 3D grid. The planning
month is therefore fixed for each path, even though paths may be longer than 1
month in duration.

The travel time from one waypoint (the source) to another (the destination),
within the layer characterising the planning month, can be found by using our
3D-DSP method to compute the travel time from the centre of the cell
containing the source, to the centre of the cell containing the destination,
and adding the travel times required to get from the source and destination to
and from their respective centre points. The entire path is constructed within
the horizontal layer. The resulting route and travel-time is stored in a
structure called the path-book, indexed by the planning month. This first
stage is called the path-book construction stage.

Given a suitably dense set of waypoints, many of the paths in the path-book
are of less than one month duration. However, the path-book construction stage
does not limit routes to be no more than one month long. Some routes can be
completed beyond the outer limits of the ice edge and therefore are not
affected by the changing sea ice extent. In cases where a long route does
visit waypoints in ice-prone areas, planning within a layer results in
unexecutable routes because the ice conditions relevant to later parts of the
routes are different from those modelled in the planning month. Even these
routes might have long sections that are independent of the ice.

The second stage runs Dijkstra’s algorithm, over a meta-level graph, to find
optimal multi-month routes by considering compositions of the routes in the
path-book. The composition process steps vertically between layers in the
3-dimensional grid, allowing different parts of the routes to be planned under
different monthly abstractions. Sometimes, several layers must be traversed in
one vertical step, corresponding to time spent waiting for a waypoint to
become accessible. We impose the constraint that waiting time cannot be spent
in ice-bound locations. This second stage is called the path-composition stage
and is described in detail in section 4.2.2.

####  4.2.1 Route-planning

We begin by describing how the horizontal routes comprising the path-book are
constructed using Dijkstra’s and Newton’s methods.

###### Definition 2

Two cells, A and B, are adjacent if the straight line connecting the centre of
A to the centre of B passes through no other cell.

A key question is how best to pass between adjacent cells in a grid defined by
a given homogeneity configuration. When all cells are the same size, a given
cell can have at most 8 adjacent cells. In the refined grid, a cell can have
as many as 36 adjacent cells when the depth of subdivisions is limited to 3 (8
on each side plus 4 corner cells). Each subdivision approximately doubles the
number of adjacency relationships that must be considered.

The path between two adjacent cells is found by optimising travel in the cell
array representation of the vector field. The non-diagonal cases to be
considered are adjacent arrangements where the cell being travelled to is at
the left, right, above or below the cell being travelled from. All these
arrangements can be normalised into an adjacent pair where the task is to go
from the left centre to the right centre. This requires rotation of the
arrangement and consequent interpretation of the $x$ and $y$ axes and the
current vectors. After normalisation, the path between the adjacent cells goes
through an intermediate point on the edge between them, chosen to optimise the
travel time between their centre points. This is of course not always the
horizontal line, because of the currents present in the two cells. We use
Newton’s method to optimise the choice of this point, which we call $yval$.

In a constant current, the fastest route between two points will be the
straight line between them. To achieve the straight line path the vehicle must
maintain a constant net velocity along that line, which will be the sum of the
current velocity and the velocity the vehicle is maintaining relative to the
water. If the vehicle crosses the central boundary between cells at $yval$
(measured relative to the central crossing point), then the time, $t$, it
takes to travel from the centre of the cell to the crossing point on the
boundary satisfies:

$t(\vec{u}+\vec{v})=\vec{d}$ (3)

where $\vec{u}$ is the current vector in the cell, $\vec{v}$ is the velocity
of the vehicle and $\vec{d}$ is the vector from the centre of the cell to the
crossing point on the boundary. The speed of the vehicle is constant on this
path and has magnitude $s=|\vec{v}|$. Rearranging Equation 3, we have:

$t\vec{v}=\vec{d}-t\vec{u}$ (4)

the square of which is:

$s^{2}t^{2}=|u|^{2}t^{2}-2Dt+|d|^{2}$ (5)

where $D=\vec{u}.\vec{d}$. This equation can be solved for $t$:

$t=\frac{\sqrt{D^{2}+|\vec{d}|^{2}(s^{2}-|\vec{u}|^{2})}-D}{s^{2}-|\vec{u}|^{2}}$
(6)

with a special case when $s^{2}=|\vec{u}|^{2}$ which is when the speed of the
current matches the speed of the vehicle. In this case,
$t=\frac{|\vec{d}|^{2}}{D}$, being undefined when $D=0$ (the case in which the
current is orthogonal to the intended direction of travel). If $D$ is
negative, this latter case is degenerate, with the vehicle not fast enough to
overcome the current in the intended direction. The solution for $t$ defines
the function $ttl(yval)$ for a given current vector, cell size and speed. The
relevant vectors are shown in Figure 4.

Figure 4: The vector diagram showing how the time taken to travel through the
left-hand cell is calculated.

Our function is then the sum of the two travel times, $ttl(yval)+ttr(yval)$ ,
which we want to minimise by choosing the best $yval$ possible. The iteration
for Newton’s method is then:

$yval_{k+1}=yval_{k}-\frac{ttl^{\prime}(yval_{k})+ttr^{\prime}(yval_{k})}{ttr^{\prime\prime}(yval_{k})+ttr^{\prime\prime}(yval_{k})}$
(7)

When the adjacent pair is diagonally arranged, the travel time between them is
initialised to be the straight line between their centres through the corner
point connecting the two cells.

In the refined parts of the grid, the fastest route between two adjacent cells
of different sizes will be the line through a $yval$ that is displaced up or
down by an amount depending on the number of subdivisions performed. The fact
that the fastest path will always be a straight line that passes through the
edge of the smaller cell is proven in Theorem 1 in Appendix A.

The details of the modifications to Equations 3\- 7, required to manage
adjacent cells of different sizes, are given in Appendix C. The method can be
called lazily to compute only the travel times required for the specific
Dijkstra computation. When the required $yvals$ have been determined,
Dijkstra’s algorithm can select the best sequence of cells to pass through en
route between two waypoints.

As well as responding to currents, planned routes should expose the vehicle to
as little risk as possible. The formal definition of risk is as the product of
the probability of an event occurring and the severity of the event. The
probability of an event occurring depends on a model of the interaction
between the vehicle and ice at different concentrations over time, and
deriving this model is beyond the scope of this paper. We therefore leave the
development of a full risk model for future work, and instead consider the
simple concepts of risk avoidance and risk exposure. Risk avoidance is
achieved during route-planning, by the route planner choosing not to enter
cells with a mean ice concentration above a given threshold, while risk
exposure is taken to be the time for which a vehicle is exposed to mean ice
concentrations above some threshold during simulation of a route. As well as
using the effect of the current to determine which adjacent cell to visit,
Dijkstra’s algorithm compares the mean ice concentration of each cell,
according to the monthly abstraction in which the route is being planned, with
the risk avoidance threshold $\phi$. If $\phi$ is exceeded, the cell is not
entered.

####  4.2.2 3D-DSP: Route composition

A vehicle moving at 3km/h in the area of study is faster than the current in
any cell through which it travels222The fastest current in our region of the
Southern ocean is the Antarctic Circumpolar Current (ACC) which, at its
fastest, travels at 40-60cms (about 2.1km/h). This allows routes to
approximate direct paths between source and destination waypoints. The effect
of the currents is seen in the selection of $yvals$ and the durations of the
journeys, but paths are not distorted by currents pushing vehicles off-course.
In contrast, vehicles travelling at 1km/h are pushed off-course by the faster
currents (such as those in the Drake Passage) if no route to avoid the
currents can be found.

Figure 5 shows a projection of a route, from Argentine Basin to Marguerite
Bay, planned using a 1km/h vehicle under the February 2013 monthly
abstraction. The marked effect of the current (shown by the arrows) can be
seen. Cells considered land are coloured green, while cells considered ice are
white or grey, depending on their ice concentration in the planning month. The
blue cells are open water.

This route was planned horizontally during the path-book construction, but is
actually of almost 6 months duration. By the time the vehicle reaches the
Northern Peninsula, the remainder of the proposed route is infeasible because
it is in ice (the red points on the route shown where ice is encountered when
the route is simulated starting in February 2013).

Figure 5: A 175-day journey planned in February 2013 for a 1km/h vehicle.

Our path composition method removes infeasible sections of routes and replaces
them with segments planned in months in which they can be completed without
running into ice. A composed route consists of segments from the path-book,
each planned in the associated month, assembled by Dijkstra’s algorithm to
achieve the shortest overall journey-time between source and destination.

The composition is performed using a graph, $G_{M}$, representing paths
available at different months in the year. $G_{M}=<V_{M},E_{M}>$ is
constructed with vertices, $V_{M}$, being the set of pairs $\langle
Waypoint,Month\rangle$, and edges, $E_{M}$, as follows. If there is a path
from wp1 to wp2 in month $m$, of at most 30 days’ duration, then $E_{M}$
contains edges from $\langle\textit{wp1},m\rangle$ to
$\langle\textit{wp2},m\rangle$ and $\langle\textit{wp1},m\rangle$ to
$\langle\textit{wp2},m+1\rangle$. If wp1 is accessible in months $m$ and $m+1$
then $E_{M}$ contains the edge $\langle\textit{wp1},m\rangle$ to
$\langle\textit{wp1},m+1\rangle$. This edge represents waiting at wp1 during
month $m$ and ensures waiting is possible only at ice-free locations.

Figure 6: The composed journey for a 1km/h vehicle, from Argentine Basin to
Marguerite Bay, leaving in December and arriving in May. The route consists of
5 subroutes, and the projection is made up from the 5 corresponding tiles.

Algorithm: Compose($init$,$goal$,$G_{M}$)

Data: A year graph, $G_{M}$, a specified source waypoint $init$, and a
specified destination waypoint $goal$.

Result: The shortest path from $init$ to $goal$ in $G_{M}$.

begin

paths $\leftarrow\emptyset$;

for  _i=Jan..Dec_ do

$\mathcal{S}\leftarrow\langle init,i\rangle$;

$\mathcal{D}\leftarrow\emptyset$;

for  _j=Jan..Dec_ do

$\mathcal{D}\leftarrow\mathcal{D}\cup\langle goal,j\rangle$;

paths $\leftarrow$ paths $\cup$ Dijkstra($G_{M},\mathcal{S},\mathcal{D}$);

return shortest path in paths;

Algorithm 2 Long paths are composed by running Dijkstra’s algorithm over the
path-book to find the shortest sequence of path components leading from $init$
to $goal$.

In the use of Dijkstra’s algorithm in Algorithm 2, we track the total days
travelled and prevent use of an edge within the same month once the arrival
day extends into the next month.

The composed version of the route shown in Figure 5 is given in Figure 6. The
ice shown in the final segment is that of the May abstraction.

While the individual sections in the path-book minimise travel-time, the
composed route may not be globally optimal as it is forced to pass between
waypoints visited in the path-book components. This can be seen in Figure 6.
Its structure is close to that of the optimal path in Figure 5 but, because it
is composed of 5 short paths, it visits a number of points en route (including
Burdwood Bank and wp8). A detailed calendar format of this path can be found
in section 5.3.

It may be necessary for the vehicle to wait en route, to allow iced-up regions
to become accessible. The value of such a journey can be increased by electing
to wait at intermediate waypoints of scientific interest or where the vehicle
is safe from environmental impacts. We provide an arbitrarily chosen
collection of anonymous waypoints (identified as wps 1-34 in Appendix B) to
serve this purpose, shown in Figure 2. We leave identification of
scientifically useful and safe way-stations for future work.

##  5 Evaluation of the Route-planning Method

Our experiments focus on waypoint accessibility, route efficiency and the
degree of exposure to risk implied by routes obtained under the different
combinations of homogeneity configuration, risk avoidance and risk exposure
settings. We contrast two speeds of vehicle: 1km/h, a speed typical of the
commercially available gliders commonly used for ocean observation, and 3km/h,
a speed more commonly associated with powered AUVs, and plausibly achievable
by next-generation gliders. In all of the results presented we focus on
1-month routes, since these are the basic elements of all of our journeys.

The number of waypoints accessible, and hence the number of journeys possible,
in each year depends on the homogeneity configuration used and how risk-averse
the route-planner is in deciding whether or not a cell is safe to enter. We
run experiments with the Weak and Strong homogeneity configurations described
in section 4.1.4. For each year, and for each combination of vehicle speed,
risk-avoidance threshold, and homogeneity configuration, we generate a path-
book consisting of twelve sets (one for each month) of optimal routes between
all pairs of waypoints.

We compare two risk avoidance thresholds for planning, $\phi=5\%$ and
$\phi=8\%$. Both are used by all three homogeneity configurations. The impact
of the threshold on waypoint accessibility is that the lower the value of
$\phi$ the more risk-averse the route-planner is and the fewer routes are
available in any year.

To evaluate the routes planned in a given configuration and year, $y_{i}$, we
simulate running the routes against the daily ice concentrations recorded in
the B-SOSE data set for a selected year $y_{j}$ ($y_{i}$ and $y_{j}$ may be
the same or different). The daily ice data is taken as the ground truth
against which we evaluate the impact of the monthly abstraction. We use a
fixed threshold of $\psi=12\%$ for determining risk exposure during
simulation.

###  5.1 Risk Exposure

To determine the extent to which the route is exposed to risk, we consider
each leg of the journey, where a leg is a constant bearing course from point
to point, either from the centre of a cell to a crossing point, or a crossing
point to a cell centre. Where a leg would start and end on different days
(based on the duration of the journey up to the start point), the leg is
subdivided to align with the day boundaries. In this way, we can associate
each leg with a single day of ice conditions. The specific subset of B-SOSE
data points that lie inside the rectangle (aligned with the latitude-longitude
grid) defined by the end points of the leg (rounded to include points around
the waypoints), as illustrated in Figure 7, is then averaged to give the
average ice conditions encountered on that leg.

Figure 7: Data points used in the evaluation of a journey leg – B-SOSE data at
sixth of a degree resolution.

As the vehicle traverses a leg in simulation, it either encounters ice above
the $\psi=12\%$ threshold or it does not. If it does, the duration of the leg
is added to the total number of hours of risk encountered so far. At the end
of the journey, the total number of risk hours is divided by 24 to give the
number of risk days encountered on the journey. The number of risk days
encountered at the $\psi=12\%$ level is taken as a measure of the risk
exposure experienced by the vehicle on the journey.

While the monthly abstraction captures how the ice extent changes between
months, it obscures detail about how the ice extent changes within months. For
example, the Weddell Sea starts January with extensive ice cover, which
quickly collapses into February. The ice extent then builds steeply in March
and April. A monthly abstraction of the sea ice in a cell into a mean value
obscures these patterns.

To explore how the ice growth and retreat processes affect exposure to risk,
we count the number of risk days that occur within the month when travel was
planned (the first 30 days of simulation). We separate our analysis into two
parts: intra-year comparisons and inter-year comparisons. The intra-year
comparisons consider routes planned under the monthly abstractions of a given
year and simulated against the daily ice of that same year. The intra-year
comparisons provide insights into the information content of the monthly
abstraction itself. The inter-year comparisons consider routes planned under
the monthly abstractions of each of the 5 years, simulated against the daily
ice encountered in each of the other years. These comparisons show how robust
the monthly abstraction is to annual variability.

Risk exposure is a consequence of the monthly abstraction, and can therefore
never be avoided in planning, because some abstraction must always be made.
Any risk exposure not accounted for by the monthly abstraction is avoided by
using the composition method when a route of duration greater than 30 days is
needed. As an optimisation, routes that are not affected by ice can be taken
straight from the path-book.

The details of the risk exposure encountered by vehicles travelling at 1km/h
and 3km/h, using the Strong homogeneity configuration and the $\phi=8$% risk-
avoidance threshold, are given in Appendix B. The plots in Figures 8 and 9
represent the results of planning routes in the identified year and simulating
them in all 5 years.

Figure 8: Risk exposure, showing routes generated in the year specified in the
label, simulated in all 5 years. Plans are constructed using the 1km/h vehicle
in Strong using $\phi=8\%$. In each plot, the blue line is the intra-year
percentage of travel days that are risk days. The black dotted line shows the
mean percentage of risk days over all other years, encountered across all
routes in each month (the inter-year comparison). The minimum and maximum
bounds on risk exposure in each month are also shown. The green line indicates
the number of travel days planned in each month.

Figure 9: Risk exposure, showing routes generated in the year specified in the
label, simulated in all 5 years. Plans are constructed using the 3km/h vehicle
in Strong using $\phi=8\%$. In each plot, the blue line is the intra-year
percentage of travel days that are risk days. The black dotted line shows the
mean percentage of risk days over all other years, encountered across all
routes in each month (the inter-year comparison). The minimum and maximum
bounds on risk exposure in each month are also shown. The green line indicates
the number of travel days planned in each month.

Figure 10: The ice extent in the Weddell Sea region in the 5 years 2013-2017.
Data from the B-SOSE reanalysis product.

It is interesting to compare the risk exposure results, shown in Appendix B,
with the ice extent in the Weddell Sea region shown in Figure 10. Considering,
for example, Figure 9, the black lines in the plots quite faithfully follow
the patterns of the corresponding years in Figure 10. In the months in which
the ice in the planning year is less extensive than in other years, the mean
risk exposure is greater, while when the ice is more extensive than in other
years the mean risk exposure is smaller. This is because, if route-planning is
done under the pessimistic assumption that the ice extent is greater than in
other years, the risk exposure encountered in simulation tends to be smaller
than if an optimistic assumption is made. As an example, consider the plot for
2016 in Figure 9, against the 2016 line in Figure 10. The 2016 line in Figure
10 tracks the minimum sea ice area for the first 90 days of the year, and the
means (and corresponding min and max bounds) seen in the 2016 plot of the
figure are high over the first 3 months. The 2016 line then tracks the maximum
ice until day 250, and the mean risk exposure in the autumn and winter months,
shown in the plot, is correspondingly low over this period. After day 250, the
2016 line descends rapidly to track the minimum ice again over the final
hundred days of the year, and the corresponding mean risk exposure is high in
the final 3 months of the plot. The ranges of the min-max bounds in each month
are largest when the gradient of the line is steepest, because the average
over a rapidly changing month is less informative. A similar pattern to that
seen in 2016 can be observed in all of the other plots of the figure. In
particular, it can be noted that 2015 is essentially the maximal ice year
according to Figure 10, and the plot for 2015 shows very low mean risk
throughout (as well as very low waypoint accessibility). The plots for 2013
and 2014 both have low mean risk exposures in the last quarter of the year,
and Figure 10 shows that they both tracked the maximum sea ice area across all
years in those months.

###  5.2 Accessibility and Efficiency

Figures 11 and 12 show how waypoint accessibility, route efficiency and risk
exposure are affected by the choice of homogeneity configuration and risk
avoidance threshold.

Figure 11: Waypoint accessibility, journey efficiency and daily risk using a
1km/h vehicle. Figure 12: Waypoint accessibility, journey efficiency and daily
risk using a 1km/h vehicle.

In both figures, the numbers of journeys and risk days are read against the
left-hand axis, and the two quotients are read against the right-hand axis.
The figures presented are totals across all 5 years. The first 3 sets of
columns relate to $\phi=5$%, the second set relates to $\phi=8$%. The number
of routes available is taken as a measure of waypoint accessibility, while the
number of travel days per journey is a measure of route efficiency. The number
of intra-year risk days encountered is used as a measure of the intrinsic
safety of the route-planning method, while the (scaled) risk/day measurement
quantifies how intrinsically risky typical days of travel are. In Figure 11,
showing the results for the 1km/h vehicle, it can be seen that there are never
more than 8000 routes available in any of the combinations of parameters
considered. The number of routes available increases with further splitting.
The average duration of a within-month journey is 18 days when a 1km/h vehicle
is used. Risk exposure is always lower when $\phi=5$% is used than when
$\phi=8$% is used, and the daily risk is also lower in the $\phi=5$%. However,
interestingly, in both risk avoidance levels, daily risk drops off slightly in
the Strong homogeneity configuration, indicating that increased splitting
allows the ice to be more effectively avoided on some routes, without loss of
accessibility or efficiency. This holds for both risk avoidance thresholds,
indicating that a strong homogeneity condition is better than a weak one and,
since there are more routes available under $\phi=8$% than under $\phi=5$%, a
less conservative risk-avoidance threshold is also better.

The same picture is seen in Figure 12, except that – predictably – the number
of routes available (and hence waypoint accessibility) is much greater when a
faster vehicle is used. The average duration of a within-month journey is 16
days for the 3km/h vehicle, but the 3km/h vehicle covers about 3 times the
distance covered by the 1km/h vehicle in the same period of time. The amount
of risk exposure encountered is slightly higher for the 3km/h vehicle, but
this is because a faster vehicle covers more ground and encounters a greater
variety of conditions.

###  5.3 Composed Routes

Routes of many months duration, covering thousands of kilometers through areas
characterised by changing seasonal ice, can be constructed using the
composition method.

Tables 1 and 2 show a selection of journeys constructed using our composition
method, using the 3km/h and 1km/h vehicles. Given our models, some of the
destination waypoints (such as Brunt and Northern Weddell Sea) can only be
entered at very limited times of year. For example, Brunt is only accessible
in February when $\phi$ is as low as 8% (a less conservative risk-avoidance
threshold would increase its availability). Even using a 3km/h vehicle, the
best journey to Brunt from Amundsen Sea ($\sim$5,000km) must start in April
and spend 7 months waiting en route to get there in time. Waiting is sometimes
unavoidable on such long journeys. The journeys shown are the shortest of all
possible paths between the specified waypoints in the given year. Figure 14,
in section 5.3, is a calendar format visualisation of the Amundsen Sea to
Brunt journey using the 3km/h vehicle. The horizontal lines in the blue column
mark the days of travel, with the visited waypoints marked on the right (their
coordinates can be found in Appendix B), and the months of travel on the left.
The yellow periods are the waiting times along the journey. Figure 13 is a
calendar format visualisation of the composed version of the journey, using a
1km/h vehicle, from Argentine Sea to Marguerite Bay, shown in Figure 5.

Route | Year | Starting | Waiting | Duration (days)
---|---|---|---|---
Amundsen Sea$\rightarrow$Brunt | 2017 | April | 7 | 328.7
Palmer$\rightarrow$South Georgia | 2015 | January | 0 | 26.8
South Sandwich Trench$\rightarrow$Elephant Island | 2014 | April | 0 | 27.5
South Georgia$\rightarrow$SR4bottom | 2017 | February | 0 | 24.5
Maud Rise$\rightarrow$Marguerite Bay | 2013 | February | 0 | 52.4
Amundsen Sea$\rightarrow$Northern Weddell Sea | 2017 | February | 0 | 63.5
Table 1: Examples of composed routes using the 3km/h vehicle. Route | Year | Starting | Waiting (months) | Duration (days)
---|---|---|---|---
Amundsen Sea$\rightarrow$Brunt | 2017 | - | - | -
Palmer$\rightarrow$South Georgia | 2015 | February | 0 | 67.4
South Sandwich Trench$\rightarrow$Elephant Island | 2014 | February | 0 | 94
South Georgia$\rightarrow$SR4bottom | 2017 | January | 0 | 82.1
Maud Rise$\rightarrow$Marguerite Bay | 2013 | February | 5 | 403
Amundsen Sea$\rightarrow$Northern Weddell Sea | 2017 | April | 6 | 346.8
Table 2: Examples of composed routes using the 1km/h vehicle.

Figure 13: Slow vehicle composed path from Argentine Basin to Marguerite Bay,
taking nearly 6 months.

Figure 14: Fast vehicle composed path from Amundsen Sea to Brunt, taking 11
months, including 6 months waiting at wp10 (waiting periods shown highlighted
in yellow).

##  6 Summary

We have developed a long-distance route-planning capability, 3D-DSP, for
underwater vehicles operating in icy regions. Our results show that:

  1. 1.

Our combination of Dijkstra’s algorithm, Newton’s method and monthly grid
refinement helps to overcome some of the limitations of existing grid-based
discretisations and wave propagation approaches without significant
computational expense.

  2. 2.

Using different combinations of risk avoidance threshold, $\phi$, and
homogeneity conditions, we can plan more or less efficient and risk-averse
journeys. Increased grid refinement can result in lower risk exposure because
of greater awareness of the locations of possible dangers.

  3. 3.

We can plan journeys of thousands of kilometers taking into account the
surface currents and the monthly behaviour of the sea ice in our segment of
the Southern Ocean.

We based this work on an historical dataset comprising 5 years of Southern
Ocean State Estimate reanalysis data. These 5 years clearly vary in terms of
sea ice extent in our region of study, as seen in Figure 10. Choosing a
maximal ice year for planning results in very little inter-year risk exposure,
but also very low waypoint accessibility. By contrast, choosing a minimal year
for planning gives high waypoint accessibility at the cost of relatively high
risk exposure. A trade-off must be made to meet operational needs.

The product of our work is a long-distance route planning capability, 3D-DSP,
analogous to that of the Tethys Long-Range AUV, for polar ocean conditions.
Our approach does not yet consider icebergs, as distinct from sea ice. This
will be addressed in future work. In addition, there are several other ideas
that remain to be explored. As mentioned, we are now developing a full risk
model to replace the risk exposure measurements. This will allow a cost
function that combines travel time with risk minimisation. Future work will
also consider other homogeneity conditions, for example: conditions based on
characteristics such as land, current, bathymetry, icebergs and other features
including the temporal variability of sea ice within months. Furthermore, as
our methods are fully generic, long-distance route-planning for different
types of vehicles, constrained by different operational limitations and
responsive to different environmental constraints, can be addressed.

## Acknowledgments

The work described in this paper has been supported by a number of grants. Fox
was supported by a Royal Society Industry Fellowship, INF\R1\201041, Brearley
by NERC Independent Research Fellowship NE/L011666/1 and also through NERC
National Capability programmes NE/R016038/1 and NE/V013254/1, and Jones by a
UKRI Future Leaders Fellowship (reference MR/T020822/1).

## Appendix A: Technical Result

The following is a result that demonstrates that the straight line path
between the centres of adjacent squares will always cross the shared edge
segment.

###### Theorem 1

Given two squares, $A$ and $B$, such that one edge of the smaller square
(without loss of generality, $B$) lies entirely within one edge of the other
square, the line connecting the centres of the two squares crosses the edge of
$B$ that is part of the edge of $A$.

Proof: Without loss of generality, take the arrangement of the squares to be
as shown in Figure 15, with the line $C_{1}C_{2}$ connecting the centres of
the two squares. The angle $\angle XC_{2}Y$ is $90\degree$. The angle $\angle
SC_{1}Y$ satisfies $0\degree\leq\angle SC_{1}Y\leq 90\degree$ because $Y$ lies
within $RS$. By complementary angles, $\angle C_{1}YX=45\degree+\angle
SC_{1}Y$ and, therefore, $45\degree\leq\angle C_{1}YX\leq 135\degree$. The
angle $\angle C_{2}YX$ is $45\degree$, so $90\degree\leq\angle C_{2}YC_{1}\leq
180\degree$. By a similar argument, $90\degree\leq\angle C_{1}XC_{2}\leq
180\degree$. Since $XY$ is contained within $RS$, the angle $\angle
XC_{1}Y\leq 90\degree$. Thus, $XC_{1}YC_{2}$ forms a convex quadrilateral and,
therefore, its diagonals intersect inside it. The intersection of the
diagonals is exactly the crossing point of the line $C_{1}C_{2}$ at the edge
of $A$ and $B$, which completes the proof.

Note that the claim remains true, trivially, when $A$ and $B$ are of equal
size and that the condition that $B$ is not degenerate (ie $X\not=Y$) ensures
that the quadrilateral $XC_{1}YC_{2}$ is also not degenerate.

Figure 15: Two adjacent squares, $A$ and $B$, with $B$ sharing its entire edge
$XY$ with edge $RS$ of $A$.

## Appendix B: Waypoints and their Coordinates

We use 29 named waypoints of scientific interest (Table 3), and an additional
34 waypoints Table 4) scattered around and between these waypoints to serve as
via points and way stations for use on long traverses.

Name | Latitude | Longitude
---|---|---
1\. Falklands | -52∘38’ 05” | -59∘52’ 48”
2\. South Georgia | -54∘52’ 45” | -37∘15’ 51”
3\. South Sandwich Trench | -56∘24’ 41” | -25∘03’ 31”
4\. Northern Peninsula | -63∘01’ 53” | -50∘56’ 24”
5\. Maud Rise | -66∘04’ 58” | -2∘49’ 39”
6\. Amundsen Sea | -69∘19’ 23” | -129∘41’ 27”
7\. Central Weddell Sea | -71∘05’ 05” | -44∘56’ 43”
8\. South Orkney Islands | -60∘29’ 10” | -45∘18’ 17”
9\. Central Scotia Sea | -58∘21’ 29” | -44∘13’ 58”
10\. Filchner Trough Overflow | -76∘59’ 34” | -35∘56’ 49”
11\. Filchner-Ronne Ice Shelf | -77∘58’ 05” | -41∘29’ 56”
12\. Argentine Basin | -43∘00’ 00” | -58∘08’ 02”
13\. Bellingshausen Sea | -70∘54’ 12” | -82∘49’ 03”
14\. Marguerite Bay | -68∘23’ 21” | -69∘14’ 37”
15\. Palmer | -64∘40’ 43” | -67∘07’ 59”
16\. Brunt | -74∘02’ 30” | -28∘17’ 05”
17\. Northern Weddell Sea | -63∘26’ 55” | -34∘12’ 07”
18\. NorthWest Georgia Rise | -52∘21’ 25” | -36∘56’ 16”
19\. Shag Rocks Passage | -52∘12’ 50” | -48∘09’ 20”
20\. Elephant Island | -60∘32’ 50” | -55∘10’ 53”
21\. Burdwood Bank | -54∘29’ 07” | -60∘08’ 20”
22\. A23 Transect (bottom) | -62∘26’ 3.07” | -33∘15’ 39.28”
23\. A23 Transect (top) | -58∘1’ 13.83” | -38∘21’ 30.47”
24\. SR4 Transect (top) | -61∘34’ 57.72” | -49∘45’ 21.68”
25\. SR4 Transect (bottom) | -68∘24’ 47.72” | -29∘18’ 47.23”
26\. Orkney Passage Entry | -60∘57’ 5.18” | -28∘34’ 43.48”
27\. Orkney Passage Exit | -60∘41’ 40.18” | -39∘10’ 34.45”
28\. Weddell Sea 1 | -70∘46’ 05” | -52∘58’ 29”
29\. Weddell Sea 2 | -66∘46’ 27” | -57∘15’ 11”
Table 3: A collection of waypoints of scientific interest located in the
region of study

Table 4 shows the 34 anonymous waypoints, and their coordinates.

Identifier | Latitude | Longitude
---|---|---
wp1 | -53∘ 45’ 0” | -32∘ 30’ 0”
wp2 | -55∘ 56’ 0” | -28∘ 7’ 0”
wp3 | -56∘ 15’ 0” | -72∘ 30’ 0”
wp4 | -56∘ 15’ 0” | -37∘ 30’ 0”
wp5 | -56∘ 18’ 0” | -55∘ 0’ 0”
wp6 | -57∘ 29’ 0” | -45∘ 0’ 0”
wp7 | -58∘ 26’ 0” | -56∘ 52’ 0”
wp8 | -58∘ 26’ 0” | -53∘ 45’ 0”
wp9 | -58∘ 26’ 0” | -53∘ 7’ 0”
wp10 | -59∘ 59’ 0” | -57∘ 22’ 0”
wp11 | -61∘ 14’ 0” | -62∘ 29’ 0”
wp12 | -61∘ 14’ 0” | -57∘ 29’ 0”
wp13 | -61∘ 14’ 0” | -52∘ 30’ 0”
wp14 | -61∘ 14’ 0” | -47∘ 30’ 0”
wp15 | -61∘ 33’ 0” | -36∘ 15’ 0”
wp16 | -62∘ 29’ 0” | -84∘ 59’ 0”
wp17 | -63∘ 26’ 0” | -37∘ 30’ 0”
wp18 | -63∘ 45’ 0” | -62∘ 29’ 0”
wp19 | -63∘ 45’ 0” | -47∘ 30’ 0”
wp20 | -63∘ 45’ 0” | -32∘ 30’ 0”
wp21 | -63∘ 45’ 0” | -22∘ 30’ 0”
wp22 | -63∘ 45’ 0” | -29∘ 59’ 0”
wp23 | -65∘ 0’ 0” | -65∘ 0’ 0”
wp24 | -65∘ 37’ 0” | -21∘ 15’ 0”
wp25 | -66∘ 15’ 0” | -42∘ 30’ 0”
wp26 | -67∘ 30’ 0” | -40∘ 0’ 0”
wp27 | -71∘ 15’ 0” | -92∘ 30’ 0”
wp28 | -71∘ 15’ 0” | -82∘ 29’ 0”
wp29 | -71∘ 15’ 0” | -77∘ 30’ 0”
wp30 | -66∘ 0’ 0” | -75∘ 0’ 0”
wp31 | -67∘ 0’ 0” | -117∘ 0’ 0”
wp32 | -67∘ 0’ 0” | -106∘ 0’ 0”
wp33 | -67∘ 0’ 0” | -102∘ 0’ 0”
wp34 | -71∘ 0’ 0” | -111∘ 0’ 0”
Table 4: A scattering of anonymous waypoints used as via points in the region
of study.

## Appendix C: Derivation of the Expression for $yval$ in the Refined Grid

Here we provide the method for deriving $yval$ when crossing between cells of
different sizes. The variable $y$ stands for the optimal crossing point value
$yval$. Figure 16 provides the context for interpreting the following
equations.

Figure 16 shows how the optimal $yval$ is shown in the case where the adjacent
cells are of different sizes. The case in which the large cell is on the left
is shown.

Figure 16: Finding the optimal $yval$ when crossing between different sized
grid cells.

Figure 16 refers to the following quantities. $t_{1}$ and $t_{2}$ are the
travel times in the left and right cells respectively. $(u_{1},v_{1})$ and
$(u_{2},v_{2})$ are the current vector components in the left and right cells
respectively. $Y$ is the vertical separation of the two centre points. $x$ is
half of the width of the left cell, and $a$ is half of the width of the right
cell. As noted above, $y$ is $yval$ itself. $b$ is the vertical distance
travelled in the right hand cell. $s$ is the speed of the vehicle.

The distances travelled in the two cells are defined as follows.

$\begin{array}[]{ll}d_{1}^{2}=x^{2}+y^{2}\\\
d_{2}^{2}=a^{2}+b^{2}=a^{2}+(Y-y)^{2}\\\ \end{array}$

The following terms are helpful in simplifying the derivations below.

$\begin{array}[]{ll}C_{1}=s^{2}-u_{1}^{2}-v_{1}^{2}\\\
C_{2}=s^{2}-u_{2}^{2}-v_{2}^{2}\\\ D_{1}=xu_{1}+yv_{1}\\\
D_{2}=au_{2}+(Y-y)v_{2}\\\ X_{1}=\sqrt{D_{1}^{2}+C_{1}d_{1}^{2}}\\\
X_{2}=\sqrt{D_{2}^{2}+C_{2}d_{2}^{2}}\\\ \end{array}$

Using these expressions, Equation 6 can be written as follows for the two
cells.

$t_{1}=\frac{X_{1}-D_{1}}{C_{1}}$ (8) $t_{2}=\frac{X_{2}-D_{2}}{C_{2}}$ (9)

Differentiation of Equation 5 with respect to $y$, in each of the two cells,
yields:

$\frac{dt_{1}}{dy}(C_{1}t_{1}+D_{1})=y-t_{1}v_{1}$ (10)
$\frac{dt_{2}}{dy}(C_{2}t_{2}+D_{2})=y-Y+t_{2}v_{2}$ (11)

The minimum travel time is achieved when
$\frac{dt_{1}}{dy}+\frac{dt_{2}}{dy}=0$. Therefore, using equations 10 and 11
we require $F(y)=0$ where:

$F(y)=X_{1}(y-Y+\frac{v_{2}(X_{2}-D_{2})}{C_{2}})+X_{2}(y-\frac{v_{1}(X_{1}-D_{1})}{C_{1}})$
(12)

Referring back to Equation 7, note that:

$F(y)=X_{1}X_{2}(ttl^{\prime}(y)+ttr^{\prime}(y))$

The derivation in section 4.2.1 is a specialisation of the above, in which
$Y=0$ and $a=x$ (the case in which the adjacent cells are the same size).

## References

  * [1] S. J. Bijlsma.  On Minimal-Time Ship Routing.  PhD thesis, University of Delft, 1975.
  * [2] K.P. Carroll, S.R. McClaran, E.L. Nelson, D.M. Barnett, D.K. Friesen, and G.N. William.  Auv path planning: an a* approach to path planning with consideration of variable vehicle speeds and multiple, overlapping, time-dependent exclusion zones.  In Proceedings of the 1992 Symposium on Autonomous Underwater Vehicle Technology, pages 79–84, 1992.
  * [3] Edsger W Dijkstra.  A note on two problems in connexion with graphs.  Numerische mathematik, 1(1):269–271, 1959.
  * [4] R. A. Finkel and J. L. Bentley.  Quad trees: a data structure for retrieval on composite keys.  Acta Informatica, 4(1):1––9, 1974.
  * [5] B. Hobson, James Bellingham, Brian Kieft, Rob McEwen, Michael Godin, and Yanwu Zhang.  Tethys-class long range auvs - extending the endurance of propeller-driven cruising auvs from days to weeks.  Proc. IEEE/OES Autonomous Underwater Vehicles, pages 1–8, 09 2012\.
  * [6] Zhiyuan Li, Jonas W. Ringsberg, and Francisco Rita.  A voyage planning tool for ships sailing between Europe and Asia via the Arctic.  Ships and Offshore Structures, 15(S1):S10–S19, 2021.
  * [7] Matthew R Mazloff, Patrick Heimbach, and Carl Wunsch.  An eddy-permitting southern ocean state estimate.  Journal of Physical Oceanography, 40:880–899, 2010.
  * [8] Opeyemi S. Orioke, Tauhidul Alam, Joseph Quinn, Ramneek Kaur, Wesam H. Alsabban, Leonardo Bobadilla, and Ryan N. Smith.  Feedback Motion Planning for Long-Range Autonomous Underwater Vehicles.  OCEANS 2019 - Marseille, OCEANS Marseille 2019, 2019-June, 2019\.
  * [9] Clément Pêtrès, Yan Pailhas, Pedro Patrón, Yvan Petillot, Jonathan Evans, and David Lane.  Path planning for autonomous underwater vehicles.  IEEE Transactions on Robotics, 23(2):331–341, 2007.
  * [10] Stephen C. Riser, Howard J. Freeland, Dean Roemmich, Susan Wijffels, Ariel Troisi, Mathieu Belbéoch, Denis Gilbert, Jianping Xu, Sylvie Pouliquen, Ann Thresher, Pierre Yves Le Traon, Guillaume Maze, Birgit Klein, M. Ravichandran, Fiona Grant, Pierre Marie Poulain, Toshio Suga, Byunghwan Lim, Andreas Sterl, Philip Sutton, Kjell Arne Mork, Pedro Joaquín Vélez-Belchí, Isabelle Ansorge, Brian King, Jon Turton, Molly Baringer, and Steven R. Jayne.  Fifteen years of ocean observations with the global Argo array.  Nature Climate Change, 6(2):145–153, 2016.
  * [11] Debabrata Sen and Chinmaya P. Padhy.  An approach for development of a ship routing algorithm for application in the north indian ocean region.  Applied Ocean Research, 50:173–191, 2015.
  * [12] Pierre Testor, Brad DeYoung, Daniel L. Rudnick, Scott Glenn, Daniel Hayes, Craig Lee, Charitha B. Pattiaratchi, Katherine L. Hill, Emma Heslop, Victor Turpin, Pekka Alenius, Carlos Barrera, John Barth, Nicholas Beaird, Guislain Becu, Anthony Bosse, François Bourrin, Alex Brearley, Yi Chao, Sue Chen, Jacopo Chiggiato, Laurent Coppola, Richard Crout, James Cummings, Beth Curry, Ruth Curry, Richard Davis, Kruti Desai, Steven DiMarco, Catherine Edwards, Sophie Fielding, Ilker Fer, Eleanor Frajka-Williams, Hezi Gildor, Gustavo Goni, Dimitri Gutierrez, Stephanie Hanson, Peter Haugan, David Hebert, Joleen Heiderich, Karen J. Heywood, Patrick Hogan, Loïc Houpert, Sik Huh, Mark E. Inall, Masso Ishii, Schin ichi Ito, Sachihiko Itoh, Sen Jan, Jan Kaiser, Johannes Karstensen, Barbara Kirkpatrick, Jody Klymak, Josh Kohut, Gerd Krahmann, Marjolaine Krug, Sam McClatchie, Frederic Marin, Elena Mauri, Avichal Mehra, Michael P. Meredith, Travis Miles, Julio Morell, Laurent Mortier, Sarah Nicholson, Joanne O’Callaghan, Diarmuid O’Conchubhair, Peter R. Oke, Enric P. Sanz, Matthew Palmer, Jong Jin Park, Leonidas Perivoliotis, Pierre Marie Poulain, Ruth Perry, Bastien Queste, Luc Rainville, Eric Rehm, Moninya Roughan, Nicholas Rome, Tetjana Ross, Simon Ruiz, Grace Saba, Amandine Schaeffer, Martha Schonau, Katrin Schroeder, Yugo Shimizu, Bernadette M. Sloyan, David Smeed, Derrick P. Snowden, Yumi Song, Sebastiaan Swart, Miguel Tenreiro, Andrew F. Thompson, Joaquin Tintore, Robert E. Todd, Cesar Toro, Hugh Venables, Stephanie Waterman, Roy Watlington, and Doug Wilson.  OceanGliders: A component of the integrated GOOS.  Frontiers in Marine Science, 6(JUL):1–32, 2019.
  * [13] A. G. Topaj, O. V. Tarovik, A. A. Bakharev, and A. A. Kondratenko.  Optimal ice routing of a ship with icebreaker assistance.  Applied Ocean Research, 86(November 2018):177–187, 2019.
  * [14] A. Verdy and M. R. Mazloff.  A data assimilating model for estimating southern ocean biogeochemistry.  Journal of Geophysical Research: Oceans, 122:6968–6988, 2017.