11institutetext: Indian Institute of Technology Kanpur, Kanpur, India.
22institutetext: 22email: {saryan, nitik<EMAIL_ADDRESS>

# Lightweight Multi-Drone Detection and 3D-Localization via YOLO††thanks: We
gratefully acknowledge the support of NVIDIA Corporation with the donation of
the Titan Xp GPU used for this research.

Aryan Sharma 11 Nitik Jain The first and second authors have contributed
equally.11 Mangal Kothari 11 0000-0003-4713-7832

###### Abstract

In this work, we present and evaluate a method to perform real-time multiple
drone detection and three-dimensional localization using state-of-the-art
tiny-YOLOv4 object detection algorithm and stereo triangulation. Our computer
vision approach eliminates the need for computationally expensive stereo
matching algorithms, thereby significantly reducing the memory footprint and
making it deployable on embedded systems. Our drone detection system is highly
modular (with support for various detection algorithms) and capable of
identifying multiple drones in a system, with real-time detection accuracy of
up to 77% with an average FPS of 332 (on Nvidia Titan Xp). We also test the
complete pipeline in AirSim environment, detecting drones at a maximum
distance of 8 meters, with a mean error of $23\%$ of the distance. We also
release the source code for the project, with pre-trained models and the
curated synthetic stereo dataset which can be found at
github.com/aryanshar/swarm-detection

###### Keywords:

Drone detection and localization Swarm detection CNN.

##  1 Introduction

Unmanned Aerial Vehicles (UAVs) have gained massive popularity in recent
years, owing to the advancements in technology and surge in the number of use-
case UAVs those include traffic management, security and surveillance, supply
of essentials, disaster management, warehouse operations etc. Drones were
initially a military, surveillance and security tool. But in the present era,
the ecosystem around UAVs has grown into a fast growing commercial and defense
markets which has in-turn drawn investments into the UAV technology, bringing
down both shape, size and the costs. Early versions of the drone were much
larger, but as time progressed, they got smaller and smarter. Consequently
with the development of small and agile drones, their applications have time
and again raised security concerns. Their increasing use in swarm systems have
also sparked another research direction in dynamic detection and localization
of multiple drones in such systems, especially for counter-drone systems.

Drone detection is essentially a subset of the widely studied object detection
problem. Though the whole paradigm of object detection has witnessed use of
various sensors (eg. RADARS, LiDARs etc) with various novel solutions, the
real breakthrough was the use of deep-learning based methodologies for object
detection and tracking. Progressively, deep learning based solutions have
improved at the task of object detection, but have also grown bulkier and have
relied heavily on bulky computing power. Thus, these existing methods found
two roadblocks en-route to being deployed on UAVs: computing hardware and
real-time inference.

Consequently, parallel to UAV technology, development of small factor
computing board and embedded computing have made it possible to deploy deep
learning models on UAVs. Onboard object detection and localization has since
been attracting traction. This survey [9] summarizes the applications and
associated challenges in video surveillance via drones, and highlights the
recent progress and issues with the whole paradigm of drone detection and
tracking. Particularly, the problem of object localization (drones in context
of this work) is computationally expensive, since it relies on extracting
features and matching them, a problem which makes the network bulkier and
inference slower.

In this work present a light weight computer vision pipeline for dynamic
detection and localization of multiple drones. We begin with briefing about
related work in this field, followed by explaining our approach. Section 4
then presents out detailed study, where we benchmark state-of-the-art object
detection method YOLO and its variants on drone dataset, thereby selecting the
most appropriate for our pipeline. Section 5 described our depth detection
methodology in detail and concludes with the comparison of estimated depth to
the ground truth. We then conclude with the summary and proposed future work
which entails extending the network for detecting long range objects.

Figure 1: Process flow of our proposed pipeline for multiple drone detection
and localization

##  2 Related Work

For the purpose of this section, drone detection is treated as a subset of
object detection and all works have been mentioned keeping in mind their
application for drone detection. Deep learning networks have increasingly been
extending the generality of object detectors. In contrast to traditional
methods in which each stage is individually hand-crafted and optimized by
classical pipelines, deep learning networks achieve superior performance by
automatically deriving each stage for feature representation and detection.

In early years, video-based object detection is by extracting discriminant
features such as Local Binary Pattern (LBP), Scale Invariant Feature Transform
(SIFT), Histogram of Oriented Gradient (HOG) and Speeded Up Robust Features
(SURF) then using these features to train the detector. Though the classical
methods work in near real-time, they were soon outperformed by their learning
based counterparts [13].

Deep learning methods dominate the current state-of-the-art methods when it
comes to object detection, but only a selection of methods qualify for real-
time applications. Initially, single stage methods such as SSD [17] were the
ones which were most qualified for real-time application, since 2-stage
methods such as faster R-CNN [22] etc. were computationally expensive. As
fewer proposal steps with hand-crafted features are involved in single-stage
methods, they are computationally less complex than multi-state approaches
that usually prioritize detection accuracy. In practice, there was active
competition between multi-stage and single-stage methods for object detection
tasks. In 2016, Redmond et al. [21] surpassed SSD in both detection speed and
accuracy with YOLOv2. A detailed survey which was referred while choosing
detection network was done by Zhao et al. [29]. CNNs based sensors have also
been used in mobile robot localization as shown in [25].

Since UAVs can support limited payload, significant efforts have been made to
develop systems which perform computation off-board and communicate in real-
time. Lee et al. [15] demonstrated a system using Faster R-CNN, moving the
computation to an off-board computing cloud, while keeping low-level object
detection and short-term navigation onboard. The research on embedded systems,
capable of deploying deep convolution nets have lead to networks created
specifically for UAVs, targeting high speed inference and low computational
costs. Mini-YOLOv3 [18] is a real-time object specifically for embedded
applications. SSD has also been deployed on drone, demonstrating fast object
detectors as shown by Budiharto et al. [6]. Autonomous detection and tracking
of high speed ground vehicle using UAV is demonstrated in [24]. Attempts have
been made to autonomously track and land a fully actuated aircraft as shown in
[4].

The work of Hassan et al. [11] extends YOLOv3 for detecting UAVs in real-time.
Deep learning based Real-Time Multiple-Object Detection and Tracking on
embedded devices has been described in the work of Hossain and Lee [12]. Novel
approach for detection and 3D localization of drones using a hybrid motion-
based object detector has been described by Srigrarom and Hoe Chew [26].
Another approach for dynamic tracking and localization of small aerial
vehicles based on computer vision has been demonstrated by Srigrarom et al.
[27].

Since vision based methods have trouble detecting multiple objects of same
kind, those tend to fail in case of UAV swarms if not tuned properly. Hence,
non-vision based method described an efficient strategy for accurate Detection
and localization of UAV Swarms [30]. Another interesting approach was
demonstrated by Carrio et al. [7], where they use depth maps instead of normal
RGB feed to detect and localize drones. To the best of authors knowledge, we
did not find work where multiple drones have been detected from single RGB
images. Further, we then extend it presenting a pipeline for depth estimation
which is then used for 3D localization.

##  3 Proposed Methodology

We describe the overall process flow before explaining the two important
modules in the future sections. Fig. 1 illustrates the steps and processes
involved in our pipeline. The stereo camera setup on UAV captures a stereo
image which is then passed on to the object detection network. The network
then outputs two images, with bounding box labels over the drone. This has
been illustrated in the green box on the left of the image, which forms our
drone detection module

The output from network is passed on to the image centroid calculator node,
which the respective centroid in image frame coordinates. In normal depth
estimation methods, features are first identified and then the respective
disparity between them gives an estimate of depth - making the overall process
slow. We use the centroid of the two images as features, saving lot of
computation. These centroids are then passed on to the stereo triangulation
node, which outputs the estimated depth using method described in Section 5.1.
This depth information is then transformed to get the $z$-coordinate of the
drone.

Having obtained the $z$-coordinate, the $x$ and $y$ coordinates are obtained
by transforming the centroid coordinates from image frame to world
coordinates. Hence, we obtain the tuple of $(x,y,z)$ coordinates which is used
for 3D localization.

##  4 Drone Detection

Given an image with a pixel grid representation, drone detection is the task
of localizing instances of drones with bounding boxes of a certain class. Fig.
2 is a graphical comparison of most of the previously mentioned models that
plots the frame rate versus the mean average precision. The frame rates
indicated are attained using a powerful Titan X GPU and far greater than we
will be able to attain, but still provide a measure for comparison among the
algorithms.

Figure 2: Comparison of different object detection algorithms tested on Pascal
VOC2007 data. Image source: https://github.com/AlexeyAB/darknet

Hence, for the purpose of drone detection, a deep convolutional neural
network-based model known as YOLO (You Only Look Once), essentially a state-
of-the-art object detection model, is chosen and trained on a dataset of drone
images. The parameters of the model have been tuned in such a way so as to
better scope down the task of object detection to that of a drone detection.
This work involves the training of four versions of the YOLO model namely,
YOLOv3 [21], YOLOv3-tiny [2], YOLOv4 [5], YOLOv4-tiny [31], and we compare
them on the basis of some performance metrics to choose the one best suited
for our problem.

###  4.1 Dataset for Training

Dataset forms an integral part of the training of a neural network as the
quality of a dataset directly impacts the model accuracy. In this work, the
dataset has been curated from two main sources: Drone-Net [8] and Mendeley
Drone Dataset [3]. In addition to this, some images have been taken from the
internet and labeled manually using Labelbox annotation-tool [14], so as to
enrich the dataset with images containing multiple drones. Other than the
images of drones, the dataset also contains images of non-drone, drone-like
“negative” objects. This is done so as to avoid our model from overfitting.

The dataset contains 5529 images along with annotated files corresponding to
each image, containing parameters of bounding box such as height, width,
center $x$, $y$ coordinates, and object class. The dataset is further divided
into training set having 5229 images and test set having 300 images. Two
different image resolutions, 450 $\times$ 280 and 1280 $\times$ 720 were
tested for model training and the resolution 1280 $\times$ 720 yielded better
results. Hence this resolution was chosen for every other version of YOLO that
was trained.

Figure 3: Sample images from Dataset used for training the detection pipeline.
Higher resolution images can be found on Github repository

###  4.2 Training

Since we were primarily comparing versions of YOLO, Darknet ([20]) has been
used to train the models. NVIDIA TitanXP (12 GiBs vRAM) was used for training
along with GPU acceleration libraries CUDA11.3 and cuDNN8.2.

Since training heavily depends on the size of dataset, transfer learning based
approach [19] has been used to counter the problem of over-fitting, which is
prominent in case of smaller datasets. We begin by initializing the pretrained
model given by Redmond [20]. The training parameters for various models are
given in Table 1.

Table 1: Parameters of CFG used for training Parameters | YOLOv3 | YOLOv3-tiny | YOLOv4 | YOLOv4-tiny
---|---|---|---|---
Width | 512 | 416 | 512 | 416
Height | 512 | 416 | 512 | 416
Batch | 64 | 64 | 64 | 64
Subdivisions | 16 | 2 | 32 | 16
Channels | 3 | 3 | 3 | 3
Momentum | 0.9 | 0.9 | 0.949 | 0.9
Decay | 3 | 3 | 3 | 3
Learning Rate | 0.0005 | 0.0005 | 0.0005 | 0.0005
max_batch | 0.001 | 0.001 | 0.001 | 0.001
Policy | steps | steps | steps | steps
Steps | 3200, 3600 | 4480, 5040 | 3200, 3600 | 4480, 5040
Scale | 0.1, 0.1 | 0.1, 0.1 | 0.1, 0.1 | 0.1, 0.1
Classes | 1 | 1 | 1 | 1

###  4.3 Results

Figure 4: Results of the trained models on test images. To the left is the
test image, with outputs from various model as one moves right. Top row
contains image with a single drone, while bottom row is for multiple drones

As mentioned in Section 4.1, we evaluated our models on the test set. Selected
results (encompassing multiple and single drone) have been shown in Fig. 4.
From the output, it can be easily seen that YOLOv4 and tiny-YOLOv4 clearly
outperform YOLOv3, since YOLOv3 fails to detect all the drones due to its
lower confidence values. Further, we evaluated the inference speeds, since
those are of paramount importance for real-time deployment. From Table 2, it
can be deduced that tiny-YOLOv4 performs the best out of every model, in terms
of average confidence and inference speed. Though, literature suggests that
YOLOv4 must outperform every other listed models here (in terms of average
confidence), but that hypothesis fails here. This is largely accredited to the
the size of the dataset that the network is trained upon. Since, YOLOv4 is a
very dense network, it also requires a vast dataset, both in terms of quality
and quantity. Although, the choice of the network we’ve made is heavily based
on the fact that we require a light weight network with low memory footprint
and high inference speeds, so we move forward with tiny-YOLOv4.

Table 2: Inference time & average confidence comparison for various models on
the test images in Fig. 4 Single Drone
---
Version | Average Confidence | Inference Time (ms)
YOLOv3 | 94% | 215.25
tiny-YOLOv3 | 90% | 196.53
YOLOv4 | 94% | 212.19
tiny-YOLOv4 | 98% | 197.27
Multiple Drones
YOLOv3 | 42% | 202.09
tiny-YOLOv3 | 59% | 190.2
YOLOv4 | 75% | 195.25
tiny-YOLOv4 | 94% | 165.07

###  4.4 Why tiny-YOLOv4 ?

The major parameters over which we have compared these four YOLO models are
the confidence of prediction, inference time, and the accuracy while detecting
multiple drones. On the basis of these results, we have observed the following
:

  * •

The newer versions of YOLO performs better overall, both in terms of accuracy
and precision.

  * •

Although, YOLOv4 is a much denser network, still tiny-YOLOv4 outperforms it
especially in terms of the confidence in multiple object detection.

  * •

YOLOv4 and tiny-YOLOv4 have much better multiple drone detection accuracy than
YOLOv3 and tiny-YOLOv3

  * •

The Prediction Time of tiny-YOLOv4 is much better than the other YOLO
versions, with also the best confidence of detection.

All of these observations motivated us to choose tiny-YOLOv4 architecture as
the baseline model for the next part of the work, i.e., Depth Estimation.

##  5 Depth Estimation

The perception of depth and its estimation from a 2D image is a very
challenging problem in the field of computer vision. Stereo vision consists of
two identical cameras placed at a baseline distance, which allows them to take
images from two distinct viewpoints. The depth in this scenario is estimated
by finding the disparity of the images of the same 3D point as demonstrated by
Acharyya et al. [1]. The work of Yamaguchi et al. [28], showed that one of the
prominent ways to evaluate depth is through stereo-matching of the local
features in the stereo-image pair and then retrieving depth through
Triangulation. In this work, we propose an alternative approach so as to cut-
down the computations involving the stereo matching process. Using tiny-YOLOv4
as the drone detection framework, we feed the stereo image frames into our
trained network to retrieve the bounding box parameters of both left-camera
and the right-camera images. Now, we approximate the output vector (bounding-
box parameters) of the YOLO network as raw local features of the target object
in the image, which makes up for the required stereo-matching and reduces the
overall computation time and intensity. The approximate evaluation of depth is
then done by using the x-coordinate of the centroid of the bounding box to
first calculate the stereo-disparity and then use the Triangulation equation
to obtain the depth.

###  5.1 Stereo Triangulation

Stereo vision adds the perception of another depth dimension to a 2D digital
image. The first step in the process is to evaluate the disparities between
the images produced by the two cameras having a focal length $f$, and which
are placed at some known baseline distance, $B$. Figure 5 shows a general
arrangement for a stereo camera setup. Here, $O_{L}$ and $O_{R}$ are the
optical centers for the left-camera and right-camera, respectively.

Figure 5: Stereo Triangulation scheme with parallel camera arrangement. Image
source: [10]

The figure also depicts the projection of a distant point object $P$ onto both
the cameras. It is clear that the center of the left-camera image clearly
differs from the center of the right-camera image. This is termed as parallax
effect and it results in relative displacement of the respective image points
when seen from different viewpoints. This displacement is termed as disparity,
$\Delta x$ which is given by:

$\Delta x=x_{L}-x_{R}$ (1)

From the property of similar triangles, the parameters like baseline $B$,
focal length $f$, disparity $\Delta x$, and the depth $Z$, can be formulated
as equality of ratios,

$\frac{Z}{B}=\frac{f}{\Delta x}$ (2)

Obtaining the depth from this equation, requires rearranging these terms,

$Z=\frac{f\times B}{\Delta x}$ (3)

It is evident from equation (3) that if the disparity term is kept constant,
then decreasing the baseline B also decreases the depth $Z$. Thus, it is
suggested that to obtain depth of objects that are far away, the baseline
should be kept large.

###  5.2 Stereo Dataset formation using AirSim

A very prominent obstacle faced throughout was the lack of dataset containing
drone images, and while also considering that physically creating a dataset by
flying the drones and capturing its images from different viewpoints is a very
time-intensive task, inspired to create the stereo-dataset in a virtual
environment, which mimics reality, by using AirSim [23] and automating the
process of creating the dataset by writing a Python script for the simulation
process. AirSim exhibits certain APIs which can be used to interact with the
vehicle in the simulation programmatically to retrieve images, get state,
control the vehicle, etc.

In order to automate the simulation process of capturing the stereo images
from different viewpoints for enriching the dataset, the python script is run
simultaneously with our custom AirSim environment. The arrangement in the
simulation environment consists of two sets of drones (as shown in Fig. 6),
one are the Target Drones (whose image is being captured) that are made to fly
up and hover at a particular height, and the other ones include a set of four
reference drones (those who capture images of target drones). These reference
drones are placed equidistant from each other, having the target drones at the
center and these are made to fly such that they capture stereo images of these
target drones from different viewpoints.

Figure 6: Custom AirSim enviroment setup used for simulation. Red-cross
denotes the Target Drone and Green-cross denotes the Reference Drones Figure
7: Sample stereo-image pair of a scene after being fed through tiny-YOLOv4
network

###  5.3 Evaluation of Approximate Depth

Equation (3) suggests that given the camera parameters, namely focal length
($f$) and baseline distance ($B$), if we are able to obtain the disparity
between the images, the parameter of depth could then be evaluated. To do so,
we feed the stereo-image pair into the trained tiny-YOLOv4 network to extract
the predicted bounding box parameters. Since the stereo-camera platform that
we have taken is a parallel camera arrangement, the disparity can be
approximated as the difference in the $x$-coordinate of the target drone’s
projection in both the left-camera and right-camera image, respectively. Now,
the $x$-coordinate of the target drone’s projection in both the cameras is
approximated as the $x$-coordinate of the centroid of the bounding-box
obtained after feeding the images into the trained tiny-YOLOv4 network (see
Figure 7). This establishes the requirement for obtaining the depth parameter.
Say, $Cx_{LeftCam}$ and $Cx_{RightCam}$ denote the $x$-coordinate of the
centroid of the bounding box parameter in the left-camera image and the right-
camera image respectively, then disparity is given as:

$\Delta x_{img}=Cx_{LeftCam}-Cx_{RightCam}$ (4)

###  5.4 Formulation for Comparison

To establish a comparison of our approximate Depth-Detection technique with
the Ground-Truth, we first generate a stereo dataset of 50 images using the
simulation script as mentioned in Section 5.2. The camera parameters used for
generating these images were, Baseline, $B$ = 1.2 $m$ and focal length, $f$ =
1.2 $m$
Since, the bounding-box parameters given as output vector by the tiny-YOLOv4
network is in terms of the pixel no., the disparity $\Delta x_{img}$ value
obtained is terms of pixel number. Thus, making appropriate transformations
using the conversion 1 $pixel$ = $2.65\times 10^{-4}$ $m$ , we get

$Z=\frac{B\times f}{\Delta x_{img}(in\;pixels)}=\frac{9070.86}{\Delta
x_{img}}(\;m)$ (5)

###  5.5 Comparison with Ground-Truth

As we have the required relation for evaluating the depth of a target drone,
we begin to compare its performance. To evaluate the performance of our depth
estimation model, the ground truth values of the target drone and the
reference drone were recorded throughout the process of simulation while
capturing the dataset of 50 stereo-pair images. Out of the 50 pairs, we
compute the values of depth for a sample size of 8 stereo-image pairs using
Eq. (5) and the results are reported in Table 3 below.

Table 3: Comparison of Estimated Z-depth and Ground Truth Depth Sample | Disparity | Z-Depth | Ground-Truth | Error
---|---|---|---|---
No. | (pixel) | (m) | (m) | %
1 | 1412 | 6.42 | 5.71 | 12.32 %
2 | 1274 | 7.12 | 5.93 | 19.97 %
3 | 1173 | 7.73 | 6.37 | 21.24 %
4 | 1104 | 8.21 | 6.54 | 25.62 %
5 | 1089 | 8.33 | 6.73 | 23.75 %
6 | 1028 | 8.82 | 6.93 | 27.25 %
7 | 898 | 10.10 | 7.59 | 32.95 %
8 | 963 | 9.42 | 7.70 | 22.27 %

It is observed from the table that the proposed depth estimation model gives
close to accurate measurement of the object depth. It is also observed that
the error in the measurement of depth from our formulated model increases with
the depth. This tells us that the model is able to deliver good results when
the target drone is in a close proximity.

##  6 Conclusion

We have presented and evaluated our approach for real-time multiple drone
detection and localization in simulation environment. The pipeline is modular
with support for various object detection algorithms depending on frame rate,
YOLO and its variants have been evaluated in this work. The modern, neural net
based tiny-YOLO v4 algorithm attains higher frame rates and detection accuracy
results than leading CPU based algorithms, and coupled with our classical
stereo triangulation based depth estimation module, can be used for 3D
localization. Further, we also release the pre-trained models, source code for
simulation and the custom stereo dataset for further developments. Since
traditional feature matching methods tend to deteriorate as images have more
negative space, our method circumvents the problem since we use centroid as a
feature for stereo triangulation. Hence as future improvements, we plan to
extend our work for achieving long range drone detection using UAV-Yolo [16].

## References

  * [1] Acharyya, A., Hudson, D., Chen, K.W., Feng, T., Kan, C.Y., Nguyen, T.: Depth estimation from focus and disparity. In: 2016 IEEE International Conference on Image Processing (ICIP). pp. 3444–3448 (2016). https://doi.org/10.1109/ICIP.2016.7532999
  * [2] Adarsh, P., Rathi, P., Kumar, M.: Yolo v3-tiny: Object detection and recognition using one stage improved model. In: 2020 6th International Conference on Advanced Computing and Communication Systems (ICACCS). pp. 687–694 (2020). https://doi.org/10.1109/ICACCS48705.2020.9074315
  * [3] Aksoy, M.a.: Drone dataset: Amateur unmanned air vehicle detection (Nov 2019), https://data.mendeley.com/datasets/zcsj2g2m4c
  * [4] Bhargavapuri, M., Shastry, A.K., Sinha, H., Sahoo, S.R., Kothari, M.: Vision-based autonomous tracking and landing of a fully-actuated rotorcraft. Control Engineering Practice 89, 113–129 (2019). https://doi.org/https://doi.org/10.1016/j.conengprac.2019.05.015, https://www.sciencedirect.com/science/article/pii/S0967066118306415
  * [5] Bochkovskiy, A., Wang, C.Y., Liao, H.Y.M.: Yolov4: Optimal speed and accuracy of object detection (2020)
  * [6] Budiharto, W., Gunawan, A.A.S., Suroso, J.S., Chowanda, A., Patrik, A., Utama, G.: Fast object detection for quadcopter drone using deep learning. In: 2018 3rd International Conference on Computer and Communication Systems (ICCCS). pp. 192–195 (2018). https://doi.org/10.1109/CCOMS.2018.8463284
  * [7] Carrio, A., Tordesillas, J., Vemprala, S., Saripalli, S., Campoy, P., How, J.P.: Onboard detection and localization of drones using depth maps. IEEE Access 8, 30480–30490 (2020). https://doi.org/10.1109/ACCESS.2020.2971938
  * [8] Chuanenlin: chuanenlin/drone-net: 2664 images of drones, labeled, with trained yolo weights. example project for my article ”tutorial: Build your custom real-time object classifier” on medium. (Aug 2018), https://github.com/chuanenlin/drone-net
  * [9] Dilshad, N., Hwang, J., Song, J., Sung, N.: Applications and challenges in video surveillance via drone: A brief survey. In: 2020 International Conference on Information and Communication Technology Convergence (ICTC). pp. 728–732 (2020). https://doi.org/10.1109/ICTC49870.2020.9289536
  * [10] Fahmy, A.: Stereo vision based depth estimation algorithm in uncalibrated rectification (2013)
  * [11] Hassan, S.A., Rahim, T., Shin, S.Y.: Real-time uav detection based on deep learning network. In: 2019 International Conference on Information and Communication Technology Convergence (ICTC). pp. 630–632 (2019). https://doi.org/10.1109/ICTC46691.2019.8939564
  * [12] Hossain, S., Lee, D.j.: Deep learning-based real-time multiple-object detection and tracking from aerial imagery via a flying robot with gpu-based embedded devices. Sensors 19(15) (2019). https://doi.org/10.3390/s19153371, https://www.mdpi.com/1424-8220/19/15/3371
  * [13] Hulens, D., Goedemé, T.: Autonomous flying cameraman with embedded person detection and tracking while applying cinematographic rules. In: 2017 14th Conference on Computer and Robot Vision (CRV). pp. 56–63. IEEE (2017)
  * [14] Labelbox: Labelbox/labelbox: Labelbox is the fastest way to annotate data to build and ship computer vision applications., https://github.com/Labelbox/labelbox
  * [15] Lee, J., Wang, J., Crandall, D., Šabanović, S., Fox, G.: Real-time, cloud-based object detection for unmanned aerial vehicles. In: 2017 First IEEE International Conference on Robotic Computing (IRC). pp. 36–43 (2017). https://doi.org/10.1109/IRC.2017.77
  * [16] Liu, M., Wang, X., Zhou, A., Fu, X., Ma, Y., Piao, C.: Uav-yolo: Small object detection on unmanned aerial vehicle perspective. Sensors 20(8) (2020). https://doi.org/10.3390/s20082238, https://www.mdpi.com/1424-8220/20/8/2238
  * [17] Liu, W., Anguelov, D., Erhan, D., Szegedy, C., Reed, S., Fu, C.Y., Berg, A.C.: Ssd: Single shot multibox detector. In: European conference on computer vision. pp. 21–37. Springer (2016)
  * [18] Mao, Q.C., Sun, H.M., Liu, Y.B., Jia, R.S.: Mini-yolov3: Real-time object detector for embedded applications. IEEE Access 7, 133529–133538 (2019). https://doi.org/10.1109/ACCESS.2019.2941547
  * [19] Pan, S.J., Yang, Q.: A survey on transfer learning. IEEE Transactions on Knowledge and Data Engineering 22(10), 1345–1359 (2010). https://doi.org/10.1109/TKDE.2009.191
  * [20] Redmon, J.: Darknet: Open source neural networks in c. http://pjreddie.com/darknet/ (2013–2016)
  * [21] Redmon, J., Farhadi, A.: Yolo9000: Better, faster, stronger (2016)
  * [22] Ren, S., He, K., Girshick, R., Sun, J.: Faster r-cnn: Towards real-time object detection with region proposal networks. Advances in neural information processing systems 28, 91–99 (2015)
  * [23] Shah, S., Dey, D., Lovett, C., Kapoor, A.: Airsim: High-fidelity visual and physical simulation for autonomous vehicles (2017)
  * [24] Shastry, A.K., Sinha, H., Kothari, M.: Autonomous detection and tracking of a high-speed ground vehicle using a quadrotor uav. In: AIAA Scitech 2019 Forum. p. 1188 (2019)
  * [25] Sinha, H., Patrikar, J., Dhekane, E.G., Pandey, G., Kothari, M.: Convolutional neural network based sensors for mobile robot relocalization. In: 2018 23rd International Conference on Methods Models in Automation Robotics (MMAR). pp. 774–779 (2018). https://doi.org/10.1109/MMAR.2018.8485921
  * [26] Srigrarom, S., Hoe Chew, K.: Hybrid motion-based object detection for detecting and tracking of small and fast moving drones. In: 2020 International Conference on Unmanned Aircraft Systems (ICUAS). pp. 615–621 (2020). https://doi.org/10.1109/ICUAS48674.2020.9213912
  * [27] Srigrarom, S., Lee, S.M., Lee, M., Shaohui, F., Ratsamee, P.: An integrated vision-based detection-tracking-estimation system for dynamic localization of small aerial vehicles. In: 2020 5th International Conference on Control and Robotics Engineering (ICCRE). pp. 152–158 (2020). https://doi.org/10.1109/ICCRE49379.2020.9096259
  * [28] Yamaguchi, K., McAllester, D., Urtasun, R.: Efficient joint segmentation, occlusion labeling, stereo and flow estimation. In: Fleet, D., Pajdla, T., Schiele, B., Tuytelaars, T. (eds.) Computer Vision – ECCV 2014. pp. 756–771. Springer International Publishing, Cham (2014)
  * [29] Zhao, Z.Q., Zheng, P., Xu, S.T., Wu, X.: Object detection with deep learning: A review. IEEE Transactions on Neural Networks and Learning Systems 30(11), 3212–3232 (2019). https://doi.org/10.1109/TNNLS.2018.2876865
  * [30] Zheng, J., Chen, R., Yang, T., Liu, X., Liu, H., Su, T., Wan, L.: An efficient strategy for accurate detection and localization of uav swarms. IEEE Internet of Things Journal pp. 1–1 (2021). https://doi.org/10.1109/JIOT.2021.3064376
  * [31] Zhu, D., Xu, G., Zhou, J., Di, E., Li, M.: Object detection in complex road scenarios: Improved yolov4-tiny algorithm. In: 2021 2nd Information Communication Technologies Conference (ICTC). pp. 75–80 (2021). https://doi.org/10.1109/ICTC51749.2021.9441643


# Shape effect on ice melting in flowing water

Rui Yang1  Christopher J. Howland1  Hao-Ran Liu2  Roberto Verzicco1,3,4
Detlef Lohse1,5<EMAIL_ADDRESS>1Physics of Fluids Group, Max Planck
Center for Complex Fluid Dynamics, and J. M. Burgers Centre for Fluid
Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands
2Department of Modern Mechanics, University of Science and Technology of
China, Hefei 230027, China 3Dipartimento di Ingegneria Industriale, University
of Rome ‘Tor Vergata’, Roma 00133, Italy
4Gran Sasso Science Institute - Viale F. Crispi, 7 67100 L’Aquila, Italy
5Max Planck Institute for Dynamics and Self-Organization, 37077 Göttingen,
Germany

###### Abstract

Iceberg melting is a critical factor for climate change, contributing to
rising sea levels and climate change. However, the shape of an iceberg is an
often neglected aspect of its melting process. Our study investigates the
influence of different ice shapes and ambient flow velocities on melt rates by
conducting direct numerical simulations. Our study focuses on the ellipsoidal
shape, with the aspect ratio as the control parameter. It plays a crucial role
in the melting process, resulting in significant variations in the melt rate
between different shapes. Without flow, the optimal shape for a minimal melt
rate is the disk (2D) or sphere (3D), due to the minimal surface area.
However, as the ambient flow velocity increases, the optimal shape changes
with the aspect ratio. We find that ice with an elliptical shape (when the
long axis is aligned with the flow direction) can melt up to 10% slower than a
circular shape when exposed to flowing water. We provide a quantitative
theoretical explanation for this optimal shape, based on the competition
between surface area effects and convective heat transfer effects. Our
findings provide insight into the interplay between phase transitions and
ambient flows, contributing to our understanding of the iceberg melting
process and highlighting the need to consider the aspect ratio effect in
estimates of iceberg melt rates.

###### keywords:

##  1 Introduction

With global warming, more icebergs are breaking off from Antarctica,
accelerating sea level rise (Scambos et al., 2017). Additionally, melting
icebergs and sea ice floes also contribute significantly to climate and
environment change, including freshwater supply (Huppert & Turner, 1978),
biological productivity (Wadham et al., 2019), and carbon sequestration
(Duprat et al., 2016), making the understanding of iceberg melting rates
crucial for comprehending the interplay between icebergs and the climate
(Cenedese & Straneo, 2023). To describe the melt rate of icebergs, different
parameterization models are proposed, such as: employing empirical relations
for turbulent heat transfer over a flat plate to predict iceberg melt rates
(Weeks & Campbell, 1973), considering meltwater-plume effect on the melt rate
with ambient flows (FitzMaurice et al., 2017), and coupling effects of
turbulence, buoyant convection, and waves (Martin & Adcroft, 2010).

Various lab experiments and direct numerical simulations have been conducted
to investigate the interaction between ice melting and ambient flow in the
laboratory scale, including ice melting in still and flowing water (Dumore et
al., 1953; Vanier & Tien, 1970; Hao & Tao, 2002; Hester et al., 2021; Weady et
al., 2022), in Rayleigh-Bénard convection (Davis et al., 1984; Dietsche &
Müller, 1985; Favier et al., 2019; Yang et al., 2023b), which involves a plate
being heated from below and a plate being cooled from above (Lohse & Xia,
2010; Ecke & Shishkina, 2023), and in vertical convection (Wang et al., 2021c;
Yang et al., 2022). The density anomaly effect in freshwater causes distinct
flow regimes and ice melting morphology under different temperatures (Wang et
al., 2021b, a, c). Salinity results in even more complex mushy layer
structures (Worster, 1997), where convection also occurs (Du et al., 2023).
These studies highlight the complex interplay between ice melting and ambient
flow, leading to distinct ice morphologies and melt rates.

Icebergs and ice floes exhibit considerable variations in shape and size
(Gherardi & Lagomarsino, 2015), with horizontal extents ranging from a few
meters to several hundred kilometers. Despite extensive research on iceberg
melt rates using models, experiments, and simulations, the potential effects
of the aspect ratio on the melt rate did not get much attention. One of the
experimental results is that the overall melt rate strongly depends on the
aspect ratio (Hester et al., 2021; Cenedese & Straneo, 2023). However, this
effect of the aspect ratio is still poorly understood. Therefore, to enhance
the accuracy of iceberg melting predictions, it is imperative to consider the
aspect ratio in models of iceberg melt rates.

In this study, we investigate the influence of the ice shape (ellipse aspect
ratio) on melt rates through numerical simulations and theoretical analysis.
We focus on the scenario of ice melting in cross-flow, as detailed in Methods,
and neglect buoyancy effects and basal melting due to our focus on the top
view of the melting process, distinguishing our work from previous studies
(Couston et al., 2021; Hester et al., 2021). Our aim is to understand how the
aspect ratio affects ice melt rates by conducting a series of numerical
simulations. Our findings reveal that the aspect ratio plays a crucial role in
the melting process, leading to substantial variations in melt rates. We
further propose a theoretical model for the dependence of melt rates on the
control parameters that agrees well with the simulation results. In this
model, two contributions to the total melt rate are distinguished - surface
contact-induced melt and advective flow-induced melt. The observed strong
aspect ratio dependence can be quantitatively understood.

##  2 Numerical method and set-up

We numerically integrate the velocity field $\boldsymbol{u}$ and the
temperature field $\theta$ according to the Navier–Stokes equations,
neglecting the effect of buoyancy induced by temperature differences. The
melting process is modeled by the phase-field method, which has been widely
used in previous studies (Favier et al., 2019; Hester et al., 2021; Couston et
al., 2021; Yang et al., 2022, 2023b). In this technique, the phase field
variable $\phi$ is integrated in time and space and smoothly transitions from
a value of $1$ in the solid to a value of $0$ in the liquid. The equations are
non-dimensionalized by inflow speed $U_{0}$ as the velocity scale, ice
effective diameter $D$ as the length scale, and the temperature difference
between the ambient flow and the ice $T_{0}$ as the temperature scale. The
non-dimensionalized quantities include three velocity components $u_{i}$ with
$i=1,2,3$, the pressure $p$, the temperature $\theta$, and the phase field
scalar $\phi$. The dimensionless governing equations read

$\displaystyle\nabla\cdot\boldsymbol{u}=0,$
$\displaystyle\partial_{t}\boldsymbol{u}+\nabla\cdot(\boldsymbol{u}\boldsymbol{u})=-\nabla
p+\frac{1}{Re}\left(\nabla^{2}\boldsymbol{u}-\frac{\phi\boldsymbol{u}}{\eta}\right),$
$\displaystyle\partial_{t}\theta+\nabla\cdot(\boldsymbol{u}\theta)=\frac{1}{{RePr}}\nabla^{2}\theta+St\partial_{t}\phi,$
$\displaystyle\frac{\partial\phi}{\partial
t}=\frac{6}{5StC{{RePr}}}\left[\nabla^{2}\phi-\frac{1}{\beta^{2}}\phi(1-\phi)(1-2\phi+C\theta)\right],$

where $\beta$ is the diffusive interface thickness, which is typically set to
be the mean grid spacing. $C$ is the phase mobility parameters related to the
Gibbs–Thompson relation. We choose $C=1$, which may overestimate the Gibbs-
Thomson effect for high curvature. Therefore we avoid too extreme values of
$\gamma$, where the high curvature regions might be inaccurate. More details
can be found in previous studies (Yang et al., 2023b; Hester et al., 2020).
Simulations are performed using the second-order staggered finite difference
code AFiD, which has been extensively validated and used to study a wide range
of turbulent flow problems (Ostilla-Mónico et al., 2015; Yang et al., 2016;
Liu et al., 2022), including phase-change problems (Yang et al., 2022, 2023b).
The phase field method is applied to model the phase-change process, which has
been widely used in previous studies (Favier et al., 2019; Hester et al.,
2021; Couston et al., 2021; Ravichandran & Wettlaufer, 2021; Yang et al.,
2023b).

Figure 1: (a). An illustration of the setup for ice melting in flowing water.
The inflow is set at the left boundary with unidirectional velocity and
uniform temperature $\theta_{0}$. (b). Zoom in on the ice object. $l_{w}$ and
$l_{l}$ represent the width and the length of the object, respectively.
(c)-(f) represent the snapshots of the temperature field of ice melting in
flowing water with different aspect ratios $\gamma$, namely $\gamma=0.25$ (c),
$0.44$ (d), $1$ (e), and $2.25$ (f). The corresponding movies are shown as
supplementary materials. (g)-(j) show the snapshots of $\partial\phi/\partial
t$ (see text for more details), corresponding to the cases in (c)-(f). The
color represents the local melt rate over the surface at different times.
Also, the corresponding plot of the shift distance $\widetilde{D}_{x}$
(normalized by corresponding $l_{l}$) of the centroids as a function of the
dimensionless time is given. The mass loss rate $\dot{m}$, normalized
$\partial\phi/\partial t$, is color-coded.

The control parameters of the system are the Reynolds number $Re$, which is
the dimensionless flow strength, the Prandtl number $Pr$, which is the ratio
between kinematic viscosity $\nu$ and thermal diffusivity $\kappa$, and the
Stefan number $St$, which is the ratio between latent heat and sensible heat,
and the aspect ratio of the initial ice shape $\gamma$, which is defined as
the ratio of its width and length and will be the focus of this article:

$\displaystyle Re$ $\displaystyle=\frac{U_{0}D}{\nu},$ $\displaystyle Pr$
$\displaystyle=\frac{\nu}{\kappa},$ $\displaystyle St$
$\displaystyle=\frac{\mathcal{L}}{c_{p}\Delta T},$ $\displaystyle\gamma$
$\displaystyle=\frac{l_{w}}{l_{l}}.$ (1)

Here, $U_{0}$ is the inflow velocity, $l_{w}$ and $l_{l}$ are the initial
width and length of the ice shape (shown in fig. 1(b)), $\mathcal{L}$ is the
latent heat, $c_{p}$ is the specific heat capacity, $\Delta T$ is the
temperature difference between inflow and the ice. Due to the large parameter
space, some of the control parameters have to be fixed in order to make the
study feasible. For the time being, we fix $Pr=7$ and $St=4$ as the values for
water at $\rm 20^{o}C$, unless specified otherwise. Later we also investigate
the effect of $Pr$ and $St$ independently. Our simulations cover a parameter
range of $0\leq Re\leq 10^{3}$ ($Re=10^{3}$ corresponds to flow speed
$U=5~{}cm/s$ and ice diameter $D=2~{}cm$) and $0.1\leq\gamma\leq 2.25$.

The flow is vertically confined by two parallel boundaries with free-slip
boundary conditions for the velocity and adiabatic for the temperature (see
fig. 1(a)). In the simulations, we prescribe an ice object with area $A_{0}$
and effective diameter $D=2\sqrt{A_{0}/\pi}$ in a domain of length $L=20D$ and
width $W_{z}=5D$. Both two-dimensional ($L/W_{z}=4$) and three-dimensional
simulations ($L/W_{z}=L/W_{x}=4$) are conducted. The long and short axis
lengths of the ellipse are $l_{l}=D/\sqrt{\gamma}$ and $l_{w}=D\sqrt{\gamma}$,
respectively. Note that the finite cross-stream dimension of the computation
domain could in principle lead to blockage effects. To ensure that this does
not affect our results, we performed a test for circle shapes at $Re=400$ with
different $W_{z}$, which produced a convergent melt rate at $W_{z}=5D$.
Initially, the simulations are run with the ice object fixed at zero velocity
and $\theta_{i}=0$, and the melting process is turned off until the flow
reaches the fully developed stage, then we turn on the melting process. The
velocity and temperature at the left and right boundary are set as uniform
$U_{0}$ and $T_{0}$ by a penalty force. The initial temperature of the ice is
set as the melting temperature $\theta_{i}=0$, and the ambient fluid is set as
the inflow temperature $\theta=1$. The position of the ice front is described
as iso-contour $\phi=1/2$ as time evolves. A grid convergence test is shown in
fig. 2.

Figure 2: Resolution convergence test. (a). The normalized area as a function
of dimensionless time for different resolutions at $Re=400$. (b). The contour
plots of the ice surface at $t=8t_{0}$ (gray dashed line in (a)). Based on
this, our final choice for $Re=400$ is $N=512$.

##  3 Shape evolution and melt rate dependence on parameters

We first investigate the evolution of ice shapes, which depend on the initial
shape and the flow rate. Temperature fields under different initial shapes for
a fixed area are shown in fig. 1(c-f) ($\gamma$ increases from (c) to (f)).
From the temperature field, one can see the complex interaction between
melting and detaching wake vortices. The ice has distinct front and back
shapes - ice at the back is more deformed than at the front, due to the
detaching vortices in the wake flow, which induces reduced and non-uniform
heat flux.

To study the local melt rate at the ice front, we further plot
$\partial\phi/\partial t$ in fig. 1(g-j), whose value is zero in the solid and
liquid phases. At the solid-liquid interface, its value provides insight into
the local melt rate. It reveals that two parts of the ice melt faster (in
reddish color): the front due to inflow, and the back due to wake flow. This
melting characteristic is similar to that of dissolution (Ristroph et al.,
2012; Mac Huang et al., 2015) and to what was formed in previous experiments
on melting spheres (Hao & Tao, 2002), where the front tends to melt faster
than the back. This occurs because the liquid at the back of the solid is
cooler than at the front due to mixing induced by the wake. We investigate the
front-to-back asymmetry in the melting by plotting the distance of the ice
centroid from its initial position over time in fig. 1(g-j). For small
$\gamma$, we observe a significant shift of the centroid to the back, while
for large $\gamma$, the larger cross-section results in stronger detaching
flow, which enhances melting at the back and thus prevents a significant
movement of the centroid (see fig. 1(f)). As melting progresses, the front of
the shape remains rounded, and different back shapes are obtained for
different $\gamma$. These shapes can be attributed to the presence of the
detaching vortices, whose shedding occurs in an oscillating manner at the top
and bottom of the body. As a result, melting is favored at the top and bottom,
creating a wedge-like shape.

Figure 3: (a). The area $A$ normalized $A_{0}$ as function of time normalized
by $t_{0}=D/U_{0}$ for different $\gamma$ and fixed $Re=10^{3}$, $Pr=7$, and
$St=4$. (b). Overall melt rate $\overline{f}/f_{0}$ (from initial to complete
melt), normalized by the overall melt rate $f_{0}$ at $\gamma=1$, as a
function of $\gamma$ for different $Re$. The inset image shows the snapshots
of temperature fields for 3D simulation results ($Re=10^{3}$). The dashed line
is the theoretical curve from eq. (7), which is expected to hold for larger
$Re$.

The complicated ice-water interaction not only affects the melting shape but
also the melt rate. Fig. 3(a) shows the normalized total area change over time
for different aspect ratio $\gamma$. For increasing $\gamma$, the melt curves
exhibit a non-monotonic trend, with the melt rate first decreasing and then
increasing. We measure the time taken for the ice to completely melt as
$t_{f}$, and take $f=1/t_{f}$ to be the mean melt rate. Fig. 3(b) shows the
normalized melt rate $f/f_{0}$ as a function of $\gamma$ for different $Re$,
where $f_{0}(Re)$ is the mean melt rate for the circles ($\gamma=1$). Here,
the same trend is observed, and $\overline{f}$ depends not only on $\gamma$
but also on $Re$. In the absence of flow ($Re=0$), the minimum melt rate
occurs at the optimal aspect ratio $\gamma_{min}=1$ for a circle shape. As
$Re$ increases, $\gamma_{min}$ decreases until it converges to around
$\gamma_{min}\approx 0.5$. The same trend is also observed for 3D simulations
of a cylinder with the cross-section of different aspect ratios.

Our findings demonstrate that ice with an elliptical shape (when the long axis
is aligned with the flow direction) can melt up to 10% slower than a circle
shape when exposed to flowing water. This previously neglected shape factor
may have implications for accurately estimating the melt rate of icebergs in
previous models that neglected it (Weeks & Campbell, 1973; Martin & Adcroft,
2010; FitzMaurice et al., 2017). The non-intuitive nature of this phenomenon
raises a further question about its physical mechanism.

##  4 Theoretical model for the melt rate

To quantitatively explain the result, we consider the total ice mass budget
for the melting process:

$\frac{dA(t)}{dt}=P(\gamma,t)v_{n},$ (2)

where $P(\gamma,t)$ is the perimeter of the ice (for 2D), and $v_{n}(t)$ is
the surface averaged melt speed. Assuming the elliptical shape is not
significantly changing with time (i.e., $\gamma(t)\approx\gamma(t=0)=\gamma$),
we have the expression for the ellipse perimeter from the elliptic integral:

$P(\gamma,t)=2\sqrt{\frac{A(t)}{\pi\gamma}}\int^{2\pi}_{0}\sqrt{1-(1-\gamma^{2})\sin^{2}\alpha}d\alpha,$
(3)

where $\alpha$ is the angle. $v_{n}(t)$ would be uniform everywhere in the
absence of flow. The presence of a flow around the body increases the
temperature gradient at the ice front which enhances the heat transfer and
also makes $v_{n}(t)$ non-uniform as $v_{n}(\alpha,t)$. The flow creates a
viscous boundary layer and a thermal boundary layer, respectively of
thicknesses $\delta_{\nu}$ and $\delta_{\theta}$. The width ratio of the
boundary layers is given by the Prandtl number; $Pr=\nu/\kappa=7$ in our case
implies $\delta_{v}>\delta_{\theta}$, as illustrated in fig. 4(a-b). In order
to understand the coupling between the shape of the object and the flow around
it, we consider the Stefan boundary condition, i.e. the (dimensionless)
surface-averaged melt speed $\tilde{v}_{n}$ is related to the surface-averaged
heat flux $\overline{Nu}$:

$\tilde{v}_{n}=\frac{v_{n}}{U_{0}}=-\frac{1}{U_{0}}\frac{\kappa
c_{p}}{\mathcal{L}}\frac{\partial T}{\partial
n}=\frac{\overline{Nu}}{StPrReA^{1/2}},$ (4)

where we can use the relation from the well-known scaling of thermal boundary
layer thickness (Meksyn, 1961; Grossmann & Lohse, 2004) for the laminar
boundary layer:

$\overline{Nu}\sim\delta_{\theta}^{-1}\sim Re_{l_{w}}^{1/2}Pr^{1/3}/C(Pr),$
(5)

with $Re_{l_{w}}=Rel_{w}/D$ the Reynolds number defined by the cross-section
length $l_{w}$, $C(Pr)$ is an infinite alternating series given in Meksyn
(1961). In the limit of large $Pr$, the series for $C(Pr)$ will converge to
$C(Pr)=1$. By substituting eq. (3) and eq. (4) into eq. (2), we obtain an ODE
which is easily solved, with the result

$A(t)=A_{0}\left(1-\frac{t}{t_{f}}\right)^{\frac{4}{3}},$ (6) $\
t_{f}^{-1}\sim P(\gamma)\gamma^{1/4}\frac{Pr^{-2/3}}{C(Pr)}St^{-1}.$ (7)

Thus the overall melt rate $\overline{f}$ depends on $\gamma$ as
$\overline{f}=t_{f}^{-1}\sim P(\gamma)\gamma^{1/4}$, where $\gamma^{1/4}$
originates from the scaling of $\overline{Nu}$, representing the advective
flow-induced melt, and $P(\gamma)$ represents the surface contact-induced
melt. The contributions of $P(\gamma)$ and $\gamma^{1/4}$ are both presented
in fig. 4(b), with $P(\gamma)$ exhibiting a symmetric curve as a function of
$\gamma$ with a minimum at $\gamma=1$, which represents melting without flow.
In contrast, $\gamma^{1/4}$ displays a monotonic increase in melt rate,
indicating that ambient flow enhances melt rate. The overall melt rate
$\overline{f}\sim P(\gamma)\gamma^{1/4}$ has a non-monotonic trend with
$\gamma$, with a minimum at $\gamma_{\rm{min,theory}}=0.48$ (which is
essentially the same as the value $\gamma_{\rm{min,sim}}\simeq 0.5$ observed
by our numerical simulations). Fig. 3(b) shows the total melt rate curve,
which agrees well with the numerical results obtained at the high $Re$ in our
numerical results. For relatively low $Re$ and even $Re=0$, the advective flow
is weak, meaning $\overline{Nu}$ is close to $1$, instead of satisfying the
boundary layer relation eq. (5).

In summary, the non-monotonic relation and shift of the minimum melt rate
point are physically explained by the competition between the melt driven by
surface contact and the melt driven by the advective flow. The former is
related to the surface area, with the minimum at $\gamma=1$, while the latter
is driven by the ambient flow, leading to a monotonic increase in melt rate
with increasing $\gamma$. We confirmed that the assumption of a fixed shape as
the ice melts is valid by comparing the measured and calculated perimeters
from eq. (3) in fig. 5(a). The ratio remains close to $1$, except for large
$\gamma$, which is attributed to the strong wake flow (see fig. 1(f)). The
deviation of the assumption at large $\gamma$ also explains the deviation in
fig. 3(b) at large $\gamma$.

Figure 4: (a). A schematic of the flow and temperature fields for different
$\gamma$. The steady outer flow consists of warm water, while the attached
flow near the body forms a boundary layer (dashed line) containing the melted
fluid. It also shows that the separation point of flow moves downstream as
$\gamma$ decreases. (b). Zoomed-in view of the velocity and temperature
boundary layers, the former defined by the velocity and the latter by the
temperature gradient. (c). The theoretical curve from eq. (7), including
$P(\gamma)$, $\epsilon^{1/4}$, and $P(\gamma)\gamma^{1/4}$ as a function of
$\gamma$. Here, Pr=7, St=4.

##  5 Melt rate scalings - comparison of theory and numerical result

In order to test the scaling law $A(t)\sim(1-t/t_{f})^{4/3}$ derived in the
previous section, we conducted simulations with varying $\gamma$ and fixed
$Re=10^{3}$, $Pr=7$, and $St=4$. The evolution of the area of the object is
depicted in fig. 5(b), along with the scaling law given by eq. (6). The trends
from simulations for different $\gamma$ all follow the $4/3$ scaling, in
agreement with the theoretical result. Our findings demonstrate that
regardless of the different physical mechanisms causing the body to ablate,
whether it be through erosion (Ristroph et al., 2012), dissolution (Mac Huang
et al., 2015), or melting, the scaling laws associated with the vanishing of
the body are analogous.

Besides the dependence of the melt rate on $\gamma$, we can also obtain the
scaling dependence of the melt rate on $Pr$ and $St$ from eq. (7). To validate
this scaling law, we performed simulations with varying $St$ in fig. 5(c)
while keeping $\gamma$ fixed at $1$ and $Pr$ fixed at $7$. At large $St$
(relevant for iceberg melting in cold water), the melt rate from simulation
results follows $\overline{f}\sim St^{-1}$, which is consistent with the
scaling law from eq. (7). At small St, the scaling deviates from $St^{-1}$
because ice melts so fast that the melted fluid is still surrounding the ice,
regardless of the ambient flow. This trend has also been observed in the study
of melting in convection (Favier et al., 2019). We also conducted simulations
with varying $Pr$ in fig. 5(c) while keeping $\gamma$ fixed as $1$ and $St$
fixed as $4$. The simulation results follow $\overline{f}\sim Pr^{-2/3}/C(Pr)$
from eq. (7), implying a simplified scaling $\overline{f}\sim Pr^{-2/3}$ for
large $Pr$, which agrees with our results. In brief, our simulation results
for varying $St$ and $Pr$ both show good agreement with the scaling law
derived in eq. (7).

Figure 5: (a). The ratio between the measured perimeter from simulations and
the ideal elliptical perimeter from eq. (3) for varying $\gamma$ at $Re=400$.
The ratio close to $1$ means that our assumption is valid. (b). The normalized
area $A/A_{0}$ as a function $1-t/t_{f}$ for varying $\gamma$ at $Re=400$. All
curves follow the $4/3$ scaling as theoretically derived in eq. (6). (c). The
melt rate as a function of $St$ for fixed $\gamma=1$ and $St=4$. The dashed
gray line represents $\overline{f}\sim St^{-1}$. (d). The melt rate as a
function of $Pr$ for fixed $\gamma=1$ and $Pr=7$. The dashed gray line
represents $\overline{f}\sim Pr^{-2/3}/C(Pr)$.

##  6 Conclusions and Outlook

Through a combination of simulations and theoretical modeling, we conducted a
comprehensive investigation of the melting dynamics of a disk/elliptical shape
of ice immersed in a cross-flow, taking into account the coupled dynamics
between flow and ice melting for various initial shapes. Our results
demonstrate that in the presence of an incoming flow, the front keeps a
rounded shape during the melting process, while the back of the ice exhibits
different shapes depending on the initial shape. Our simulations reveal that
the presence of detaching vortices behind the ice plays a crucial role in
shaping the melting pattern of the object.

Furthermore, we observed that the shape of the ice also has a significant
effect on the melt rate. In the absence of flow, the circular shape
($\gamma=1$) melts more slowly than other elliptical shapes $\gamma\neq 1$.
However, in the presence of an external flow, some elliptic cylindrical bodies
($\gamma<1$) melt less rapidly than circular cylinders with the same volume.
This optimal aspect ratio depends on the flow strength, represented by $Re$.
Our physical understanding of this dependence of the melt rate on the initial
aspect ratio comes from the competition between the surface-contact-induced
melt ($\sim P(\gamma)$) and the advective-flow-induced melt
($\sim\gamma^{1/4}$). By assuming a laminar boundary layer, we derived a model
(eq. (6)and eq. (7)), accurately predicting the overall melt rate and its
scaling as a function of $\gamma$, $Pr$, and $St$. Our findings provide
insight into the rich coupling dynamics between the ice-water interface and
ambient flows and demonstrate the importance of considering ice shape in
predicting melting rates.

The approach employed in this study, which combines fully resolved direct
numerical simulations and theoretical explanations, offers the possibility to
explain other phase-change problems coupled with advective flows. We caution
that our findings reveal only a subset of the numerous factors that influence
the ice–water system dynamics. In future investigations, it would be
worthwhile to further explore additional factors such as the buoyancy force
(Favier et al., 2019; Couston et al., 2021; Wang et al., 2021b; Yang et al.,
2023b), melting near sidewalls or basal walls in the presence of ambient flow,
and the effect of dissolved salt (Huppert & Turner, 1978; Yang et al., 2023a;
Du et al., 2023). These topics have significant relevance to the modeling of
geophysical and climatological large-scale processes. Additionally, we plan to
investigate the interaction between the motion of solid objects and phase
change, as this holds significant potential for further exploration since
icebergs can move and rotate.

## Funding

We acknowledge PRACE for awarding us access to MareNostrum in Spain at the
Barcelona Computing Center (BSC) under the project 2020235589 and project
2021250115 and the Netherlands Center for Multiscale Catalytic Energy
Conversion (MCEC). We also acknowledge the support by the Priority Programme
SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaf. This
research was supported in part by the National Science Foundation under Grant
No. NSF PHY-1748958.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Cenedese & Straneo (2023) Cenedese, C. & Straneo, F. 2023 Icebergs Melting. Annu. Rev. Fluid Mech. 55, 377–402.
  * Couston et al. (2021) Couston, L.-A., Hester, E., Favier, B., Taylor, J. R., Holland, P. R. & Jenkins, A. 2021 Topography generation by melting and freezing in a turbulent shear flow. J. Fluid Mech. 911, A44.
  * Davis et al. (1984) Davis, S. H., Müller, U. & Dietsche, C. 1984 Pattern selection in single-component systems coupling bénard convection and solidification. J. Fluid Mech. 144, 133–151.
  * Dietsche & Müller (1985) Dietsche, C. & Müller, U. 1985 Influence of bénard convection on solid–liquid interfaces. J. Fluid Mech. 161, 249–268.
  * Du et al. (2023) Du, Y., Wang, Z., Jiang, L., Calzavarini, E. & Sun, C. 2023 Sea water freezing modes in a natural convection system. J. Fluid Mech. 960, A35.
  * Dumore et al. (1953) Dumore, J., Merk, H. & Prins, J. 1953 Heat transfer from water to ice by thermal convection. Nature 172 (4375), 460–461.
  * Duprat et al. (2016) Duprat, L. P., Bigg, G. R. & Wilton, D. J. 2016 Enhanced Southern Ocean marine productivity due to fertilization by giant icebergs. Nat. Geosci. 9 (3), 219–221.
  * Ecke & Shishkina (2023) Ecke, R. E. & Shishkina, O. 2023 Turbulent rotating Rayleigh–Bénard convection. Annu. Rev. Fluid Mech. 55, 603–638.
  * Favier et al. (2019) Favier, B., Purseed, J. & Duchemin, L. 2019 Rayleigh–Bénard convection with a melting boundary. J. Fluid Mech. 858, 437–473.
  * FitzMaurice et al. (2017) FitzMaurice, A., Cenedese, C. & Straneo, F. 2017 Nonlinear response of iceberg side melting to ocean currents. Geophys. Res. Lett 44 (11), 5637–5644.
  * Gherardi & Lagomarsino (2015) Gherardi, M. & Lagomarsino, M. C. 2015 Characterizing the size and shape of sea ice floes. Sci. Rep. 5 (1), 1–11.
  * Grossmann & Lohse (2004) Grossmann, S. & Lohse, D. 2004 Fluctuations in turbulent Rayleigh–Bénard convection: the role of plumes. Phys. Fluids 16 (12), 4462–4472.
  * Hao & Tao (2002) Hao, Y. L. & Tao, Y. X. 2002 Heat transfer characteristics of melting ice spheres under forced and mixed convection. J. Heat Transfer 124 (5), 891–903.
  * Hester et al. (2020) Hester, E. W., Couston, L.-A., Favier, B., Burns, K. J. & Vasil, G. M. 2020 Improved phase-field models of melting and dissolution in multi-component flows. Proc. R. Soc. A 476 (2242), 20200508.
  * Hester et al. (2021) Hester, E. W., McConnochie, C. D., Cenedese, C., Couston, L.-A. & Vasil, G. 2021 Aspect ratio affects iceberg melting. Phys. Rev. Fluid 6 (2), 023802.
  * Huppert & Turner (1978) Huppert, H. E. & Turner, J. S. 1978 On melting icebergs. Nature 271 (5640), 46–48.
  * Liu et al. (2022) Liu, H.-R., Chong, K. L., Yang, R., Verzicco, R. & Lohse, D. 2022 Heat transfer in turbulent Rayleigh–Bénard convection through two immiscible fluid layers. J. Fluid Mech. 938, A31.
  * Lohse & Xia (2010) Lohse, D. & Xia, K. 2010 Small-scale properties of turbulent Rayleigh-Bénard convection. Annu. Rev. Fluid Mech. 42 (1), 335–364.
  * Mac Huang et al. (2015) Mac Huang, J., Moore, M. N. J. & Ristroph, L. 2015 Shape dynamics and scaling laws for a body dissolving in fluid flow. J. Fluid Mech. 765, R3.
  * Martin & Adcroft (2010) Martin, T. & Adcroft, A. 2010 Parameterizing the fresh-water flux from land ice to ocean with interactive icebergs in a coupled climate model. Ocean Model. 34 (3-4), 111–124.
  * Meksyn (1961) Meksyn, D. 1961 New methods in laminar boundary-layer theory. Pergamon.
  * Ostilla-Mónico et al. (2015) Ostilla-Mónico, R., Yang, Y., Van Der Poel, E. P., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for direct numerical simulation of scalar turbulence. J. Comput. Phys. 301, 308–321.
  * Ravichandran et al. (2022) Ravichandran, S., Toppaladoddi, S. & Wettlaufer, J. S. 2022 The combined effects of buoyancy, rotation, and shear on phase boundary evolution. J. Fluid Mech. 941, A39.
  * Ravichandran & Wettlaufer (2021) Ravichandran, S. & Wettlaufer, J. S. 2021 Melting driven by rotating Rayleigh–Bénard convection. J. Fluid Mech. 916, A28.
  * Ristroph et al. (2012) Ristroph, L., Moore, M. N., Childress, S., Shelley, M. J. & Zhang, J. 2012 Sculpting of an erodible body by flowing water. Proc. Natl. Acad. Sci. 109 (48), 19606–19609.
  * Scambos et al. (2017) Scambos, T. A., Bell, R. E., Alley, R. B., Anandakrishnan, S., Bromwich, D., Brunt, K., Christianson, K., Creyts, T., Das, S., DeConto, R. & others 2017 How much, how fast?: A science review and outlook for research on the instability of Antarctica’s Thwaites Glacier in the 21st century. Glob. Planet. Change 153, 16–34.
  * Vanier & Tien (1970) Vanier, C. R. & Tien, C. 1970 Free convection melting of ice spheres. AIChE J. 16 (1), 76–82.
  * Wadham et al. (2019) Wadham, J. L., Hawkings, J., Tarasov, L., Gregoire, L., Spencer, R., Gutjahr, M., Ridgwell, A. & Kohfeld, K. 2019 Ice sheets matter for the global carbon cycle. Nat. Comm. 10 (1), 3567.
  * Wang et al. (2021a) Wang, Z., Calzavarini, E. & Sun, C. 2021a Equilibrium states of the ice-water front in a differentially heated rectangular cell (a). Europhys. Lett. 135 (5), 54001.
  * Wang et al. (2021b) Wang, Z., Calzavarini, E., Sun, C. & Toschi, F. 2021b How the growth of ice depends on the fluid dynamics underneath. Proc. Natl. Acad. Sci. 118, 10\.
  * Wang et al. (2021c) Wang, Z., Jiang, L., Du, Y., Sun, C. & Calzavarini, E. 2021c Ice front shaping by upward convective current. Phys. Rev. Fluid 6 (9), L091501.
  * Weady et al. (2022) Weady, S., Tong, J., Zidovska, A. & Ristroph, L. 2022 Anomalous convective flows carve pinnacles and scallops in melting ice. Phys. Rev. Lett. 128 (4), 044502.
  * Weeks & Campbell (1973) Weeks, W. F. & Campbell, W. J. 1973 Icebergs as a fresh-water source: an appraisal. J. Glaciol. 12 (65), 207–233.
  * Worster (1997) Worster, M. G. 1997 Convection in mushy layers. Annu. Rev. Fluid Mech. 29 (1), 91–122.
  * Yang et al. (2022) Yang, R., Chong, K. L., Liu, H.-R., Verzicco, R. & Lohse, D. 2022 Abrupt transition from slow to fast melting of ice. Phys. Rev. Fluid 7 (8), 083503.
  * Yang et al. (2023a) Yang, R., Howland, C. J., Liu, H.-R., Verzicco, R. & Lohse, D. 2023a Ice melting in salty water: layering and non-monotonic dependence on the mean salinity. J. Fluid Mech. 969, R2.
  * Yang et al. (2023b) Yang, R., Howland, C. J., Liu, H.-R., Verzicco, R. & Lohse, D. 2023b Morphology evolution of a melting solid layer above its melt heated from below. J. Fluid Mech. 956, A23.
  * Yang et al. (2016) Yang, Y., Verzicco, R. & Lohse, D. 2016 From convection rolls to finger convection in double-diffusive turbulence. Proc. Nat. Acad. Sci. 113, 69–73.