# Efficient Parallelization of a Ubiquitous Sequential Computation

Franz A. Heinsen
<EMAIL_ADDRESS>

(October 27, 2023)

###### Abstract

We find a succinct expression for computing the sequence
$x_{t}=a_{t}x_{t-1}+b_{t}$ in parallel with two prefix sums, given
$t=(1,2,\dots,n)$, $a_{t}\in\mathbb{R}^{n}$, $b_{t}\in\mathbb{R}^{n}$, and
initial value $x_{0}\in\mathbb{R}$. On $n$ parallel processors, the
computation of $n$ elements incurs $\mathcal{O}(\log n)$ time and
$\mathcal{O}(n)$ space. Sequences of this form are ubiquitous in science and
engineering, making efficient parallelization useful for a vast number of
applications. We implement our expression in software, test it on parallel
hardware, and verify that it executes faster than sequential computation by a
factor of $\frac{n}{\log n}$.111Source code for replicating our results is
available online at https://github.com/glassroom/heinsen_sequence.

##  1 Summary

Sequences of the form $x_{t}=a_{t}x_{t-1}+b_{t}$ are ubiquitous in science and
engineering. For example, in the natural sciences, such sequences can model
quantities or populations that decay or grow by a varying rate $a_{t}>0$
between net inflows or outflows $b_{t}$ at each time step $t$. In economics,
such sequences can model investments that earn a different rate of return
$a_{t}=(1+r_{t})$ between net deposits or withdrawals $b_{t}$ over each time
period $t$. In engineering applications, such sequences are often low-level
components of larger models, e.g., linearized recurrent neural networks whose
layers decay token features in a sequence of tokens.

Given a finite sequence $x_{t}=a_{t}x_{t-1}+b_{t}$ with $n$ steps,
$t=(1,2,\dots,n)$, where $a_{t}\in\mathbb{R}^{n}$, $b_{t}\in\mathbb{R}^{n}$,
and initial value $x_{0}\in\mathbb{R}$, it’s not immediately obvious how one
would compute all elements in parallel, because each element is a non-
associative transformation of the previous one. In practice, we routinely see
software code that computes sequences of this form one element at a time.

The vector $\log x_{t}$ is computable as a composition of two cumulative, or
prefix, sums, each of which is parallelizable:

$\log x_{t}=a^{*}_{t}+\log\left(x_{0}+b^{*}_{t}\right)\\\ $ (1)

where $a^{*}_{t}$ and $b^{*}_{t}$ are the two prefix sums:

$\displaystyle a^{*}_{t}$
$\displaystyle=\overset{\text{{cum}}}{\sum\limits_{t}}\log a_{t}$ (2)
$\displaystyle b^{*}_{t}$
$\displaystyle=\overset{\text{{cum}}}{\sum\limits_{t}}e^{\log
b_{t}-a^{*}_{t}}.$

The operator $\overset{\text{{cum}}}{\sum}$ computes a vector whose elements
are a prefix sum, i.e., a cumulative sum.

We obtain $x_{t}$ with elementwise exponentiation:

$x_{t}=e^{a^{*}_{t}+\log(x_{0}+b^{*}_{t})}.\\\ $ (3)

Prefix sums are associative,222 Given a sequence $a,b,c$,
$\overset{\text{{cum}}}{\sum}\Big{(}a,\overset{\text{{cum}}}{\sum}\Big{(}b,c\Big{)}\Big{)}=\overset{\text{{cum}}}{\sum}\Big{(}\overset{\text{{cum}}}{\sum}\Big{(}a,b\Big{)},c\Big{)}.$
making it possible to compute them by parts in parallel. Well-known parallel
algorithms for efficiently computing the prefix sum of a sequence with $n$
elements incur $\mathcal{O}(\log n)$ time and $\mathcal{O}(n)$ space on $n$
parallel processors Ladner and Fischer (1980) Hillis and Steele (1986). Prefix
sums generalize to any binary operation that is associative, making them a
useful primitive for many applications Blelloch (1990a) and data-parallel
models of computation Blelloch (1990b). Many software frameworks for numerical
computing provide efficient parallel implementations of the prefix sum.

The computation of two prefix sums has the same computational complexity on
$n$ parallel processors as a single prefix sum: $\mathcal{O}(\log n)$ time and
$\mathcal{O}(n)$ space. The computation of $n$ elementwise operations (e.g.,
logarithms and exponentials) on $n$ parallel processors incurs constant time
and, if the computation is in situ, no additional space.

If any $a_{t}<0$, any $b_{t}<0$, or $x_{0}<0$, one or more of the logarithms
computed in the interim will be in $\mathbb{C}$, but all elements of $x_{t}$
will always be in $\mathbb{R}$, because they are defined as multiplications
and additions of previous elements in $\mathbb{R}$, which is closed under both
operations.

##  2 Compared to Blelloch’s Formulation

Blelloch’s formulation for computing first-order linear recurrences as a
composition of prefix sums Blelloch (1990a) is more general, expressed in
terms of a binary operator $\oplus$ that is associative and a second binary
operator $\otimes$ that either is associative or can be transformed into an
associative one via the application of a third binary operator.

Our formulation applies only to the most common case, real numbers, with
scalar sum and multiplication as the first and second operators, making each
step non-associative. We find a succinct, numerically stable expression that
is readily implementable with widely available, highly-optimized
implementations of the prefix sum.

##  3 Proof

We are computing $x_{t}=a_{t}x_{t-1}+b_{t}$, for $t=(1,2,\dots,n)$, with
$a_{t}\in\mathbb{R}^{n}$, $b_{t}\in\mathbb{R}^{n}$, and initial value
$x_{0}\in\mathbb{R}$. Expand the expression that computes each element,
$x_{1},x_{2},\dots,x_{n}$, to make it a function of $x_{0}$ and all trailing
elements of $a_{t}$ and $b_{t}$, and factor out all trailing coefficients:

$\displaystyle x_{1}$ $\displaystyle=a_{1}x_{0}+b_{1}$ (4)
$\displaystyle=a_{1}\left(x_{0}+\frac{b_{1}}{a_{1}}\right)$ $\displaystyle
x_{2}$ $\displaystyle=a_{2}x_{1}+b_{2}$
$\displaystyle=a_{1}a_{2}\left(x_{0}+\frac{b_{1}}{a_{1}}+\frac{b_{2}}{a_{1}a_{2}}\right)$
$\displaystyle\vdots$ $\displaystyle x_{n}$ $\displaystyle=a_{n}x_{n-1}+b_{n}$
$\displaystyle=\Big{(}\\!\prod\limits_{t}a_{t}\\!\Big{)}\bigg{(}\\!x_{0}+\frac{b_{1}}{a_{1}}+\frac{b_{2}}{a_{1}a_{2}}+\dots+\frac{b_{n}}{{\scriptstyle\prod\limits_{t}}a_{t}}\\!\bigg{)}.$

Combine all expressions in (4) into one expression that computes all elements
of vector $x_{t}$:

$\displaystyle x_{t}$
$\displaystyle=\bigg{(}\overset{\text{{cum}}}{\prod\limits_{t}}a_{t}\bigg{)}\odot\Bigg{(}x_{0}+\overset{\text{{cum}}}{\sum\limits_{t}}\frac{b_{t}}{\overset{\text{{cum}}}{\scriptstyle\prod\limits_{t}}a_{t}}\Bigg{)}$
(5)
$\displaystyle=\bigg{(}\overset{\text{{cum}}}{\prod\limits_{t}}a_{t}\bigg{)}\odot\Bigg{(}x_{0}+\overset{\text{{cum}}}{\sum\limits_{t}}\exp\Bigg{(}\log\frac{b_{t}}{\overset{\text{{cum}}}{\scriptstyle\prod\limits_{t}}a_{t}}\Bigg{)}\\!\\!\Bigg{)}$
$\displaystyle=\bigg{(}\overset{\text{{cum}}}{\prod\limits_{t}}a_{t}\bigg{)}\odot\Bigg{(}x_{0}+\overset{\text{{cum}}}{\sum\limits_{t}}e^{\log
b_{t}-\overset{\text{{cum}}}{\sum}_{t}\log a_{t}}\Bigg{)},$

where the operators $\overset{\text{{cum}}}{\prod}$ and
$\overset{\text{{cum}}}{\sum}$ compute vectors whose elements are,
respectively, a cumulative product and sum, and $\odot$ denotes an elementwise
or Hadamard product.

Taking the logarithm on both sides, we obtain:

$\log x_{t}=\underbrace{\overset{\text{{cum}}}{\sum\limits_{t}}\log
a_{t}}_{\textstyle
a^{*}_{t}}+\log\Bigg{(}x_{0}+\underbrace{\overset{\text{{cum}}}{\sum\limits_{t}}e^{\log
b_{t}-\overbrace{\scriptstyle\overset{\text{{cum}}}{\sum}_{t}\log
a_{t}}^{\textstyle a^{*}_{t}}}}_{\textstyle b^{*}_{t}}\Bigg{)},\\\ $ (6)

which is the same as (1).

##  4 Implementation

Figure 1: Time to compute $n$ elements sequentially, relative to parallel
computation, on an Nvidia GPU. Each point is the mean of 30 runs.

We implement (3) in software. For numerical stability and slightly improved
efficiency, we modify the computation of $x_{t}$ as follows:

$x_{t}=e^{a^{*}_{t}+\operatorname{tail}(\operatorname{LCSE}(\operatorname{cat}(\log
x_{0},\log b_{t}-a^{*}_{t})))},$ (7)

where $\operatorname{cat}(\cdot)$ denotes concatenation,
$\operatorname{tail}(\cdot)$ removes its argument’s first element, and

$\operatorname{LCSE}(\cdot):=\log\overset{\text{{cum}}}{\sum}\exp(\cdot),$ (8)

commonly provided as the “LogCumSumExp” function by software frameworks for
numerical computing, applying the familiar log-sum-exp trick as necessary for
numerical stability, and delegating parallel computation of the internal
prefix sum to a highly-optimized implementation.

We test our implementation on parallel hardware and verify that it executes
faster than sequential computation by a factor of $\frac{n}{\log n}$ (Figure
1).

## References

  * Blelloch (1990a) Guy E. Blelloch. 1990a.  Prefix sums and their applications.  Technical Report CMU-CS-90-190, School of Computer Science, Carnegie Mellon University.
  * Blelloch (1990b) Guy E. Blelloch. 1990b.  Vector Models for Data-Parallel Computing.  MIT Press, Cambridge, MA.
  * Hillis and Steele (1986) W. Daniel Hillis and Guy L. Steele. 1986.  Data parallel algorithms.  Commun. ACM 29(12):1170–1183.
  * Ladner and Fischer (1980) Richard E. Ladner and Michael J. Fischer. 1980.  Parallel prefix computation.  J. ACM 27(4):831–838.


# Uplink Achievable Rate of Intelligent Reflecting Surface-Aided Millimeter-
Wave Communications with Low-Resolution ADC and Phase Noise

Kangda Zhi, Cunhua Pan, Hong Ren and Kezhi Wang (Corresponding author: Cunhua
Pan). K. Zhi, C. Pan are with the School of Electronic Engineering and
Computer Science at Queen Mary University of London, London E1 4NS, U.K.
(e-mail: k.zhi, c.pan@qmul.ac.uk). H. Ren is with the National Mobile
Communications Research Laboratory, Southeast University, Nanjing 210096,
China. (hren@seu.edu.cn). K. Wang is with Department of Computer and
Information Sciences, Northumbria University, UK. (e-mail:
kezhi.wang@northumbria.ac.uk).

###### Abstract

In this paper, we derive the uplink achievable rate expression of intelligent
reflecting surface (IRS)-aided millimeter-wave (mmWave) systems, taking into
account the phase noise at IRS and the quantization error at base stations
(BSs). We show that the performance is limited only by the resolution of
analog-digital converters (ADCs) at BSs when the number of IRS reflectors
grows without bound. On the other hand, if BSs have ideal ADCs, the
performance loss caused by IRS phase noise is constant. Finally, our results
validate the feasibility of using low-precision hardware at the IRS when BSs
are equipped with low-resolution ADCs.

###### Index Terms:

Intelligent reflecting surface, reconfigurable intelligent surface, low-
resolution ADC, hardware impairment, phase noise

##  I Introduction

Millimeter-wave (mmWave) technology will play an important role in future
networks, because of its capability of achieving high system capacity,
increased security and reduced interference[1]. However, blockage issue needs
to be tackled before the commercial application of this technology. In
specific, mmWave frequencies are susceptible to the blockage, which means
mmWave communications are difficult to be applied in urban areas with dense
buildings and vehicles.

To cope with this issue, intelligent reflecting surface (IRS), also referred
to as reconfigurable intelligent surface (RIS), has been proposed as an
attractive technology to create virtual line-of-sight (LoS) links for mmWave
communication systems[2]. IRS is composed of a large number of reconfigurable
passive elements, which is able to achieve passive beamforming by adjusting
the phase of impinging signal. Moreover, since IRS can avoid the use of radio-
frequency (RF) chains, it can effectively reduce the hardware cost and energy
consumption compared with conventional active relay nodes [3, 4, 5, 6].

Since realizing fully digital architecture in mmWave systems requires
prohibitive hardware cost and power consumption[1], IRS-aided mmWave systems
with the hybrid analog-digital architecture have been investigated in some
prior works[7, 8, 9]. Specifically, the authors in [7] demonstrated the
promising reliability and connectivity of IRS-aided mmWave systems with random
blockages. The joint passive and active beamforming optimization for single-
user mmWave networks was studied in [8], in which closed-form solution was
obtained for the single-IRS case and low-complexity iterative algorithm was
proposed for the multiple-IRS case. The authors in [9] jointly optimized
hybrid precoding at BSs and phase shifting at the IRS to minimize the mean-
squared error (MSE) for multiuser mmWave systems.

In addition to the hybrid analog-digital precoding, another promising solution
for reducing RF complexity in mmWave systems is low-resolution digital-to-
analog converters (DACs) and analog-to-digital converters (ADCs)[1]. Moreover,
compared with the DACs at the transmitter, much more power consumption is
needed for the ADCs at the receiver side [10]. Therefore, it is meaningful to
study the performance of IRS-aided mmWave systems with low-resolution ADCs.
However, to the best of our knowledge, none of the work has considered this
scenario. On the other hand, IRS phase noise, caused by its discrete phase
shifts, has been recognized as an essential factor in practical IRS-aided
systems[11, 12, 13, 14, 15]. Considering the phase estimation and quantization
errors, the authors in [11] derived the distribution of signal-to-noise ratio
(SNR) and the average error probability in single-antenna single-user systems.
In downlink multiple-input single-output (MISO) systems, the effects of phase
noise at IRS and RF impairments at BSs were studied in [12]. Under the
assumption of Rician fading channel model, the minimum number of quantization
phase bits to ensure a performance loss threshold was derived in [13] by using
the upper bound of the ergodic spectral efficiency. The impact of phase error
on the rate loss was quantified in [14] for an IRS-aided single-antenna system
with the existence of direct link. Considering the phase errors at the IRS and
the distortion noise at the BS, the authors in [15] investigated the rate
performance of IRS-aided networks and then theoretically compared it with the
relay-aided networks.

Motivated by the above discussions, we focus on the uplink IRS-aided mmWave
systems with the existence of both low-resolution ADCs and IRS phase noise,
where the IRS is used to create an additional virtual LoS link for the blocked
user. We aim to theoretically evaluate the capacity performance of this system
and characterize the joint impacts of low-resolution ADCs and IRS phase noise.
Specifically, we first derive the closed-form expression of the uplink
achievable rate. Then, based on the derived expression, we analyze the
performance loss caused by the phase noise at the IRS and quantization error
at the BS. Our asymptotic results reveal that when the number of IRS phase
shifts approaches infinity, the achievable rate is limited by the precision of
the ADCs at BSs. In addition, our results demonstrate the feasibility of
utilizing simple IRSs with a few phase shift levels when BSs are equipped with
low-resolution ADCs.

##  II System Model

Figure 1: An IRS-assisted uplink mmWave communication system.

We consider the uplink mmWave communication of an IRS-aided system with an
$M$-antenna BS, an IRS consisting of $N$ reflecting elements and a single-
antenna user, as shown in Fig. 1. Considering the blockage sensitivity in the
mmWave frequencies, we assume that the direct link is blocked by the building.
To improve the quality of service, an IRS is deployed in a proper place such
that there exists an alternative user-IRS-BS link.

We assume that both the IRS and the BS are equipped with the uniform square
planar array (USPA), and the BS can control the phase shift of each IRS’s
reflecting element. The phase shift matrix of the IRS can then be given by

$\displaystyle{\bf\Theta}={\rm{diag}}\left\\{{{\zeta_{1}}{e^{j{\theta_{1}}}},{\zeta_{2}}{e^{j{\theta_{2}}}},...,{\zeta_{N}}{e^{j{\theta_{N}}}}}\right\\},$
(1)

where ${{\zeta_{n}}}\in\left[{0,1}\right]$ and
$\theta_{n}\in\left[{0,2\pi}\right),\forall n$ denote the reflection amplitude
and the phase shift induced by the $n$th reflecting element of the IRS,
respectively. As the IRS is designed to enhance the reflected signal’s
quality, we set ${{\zeta_{n}}}=1,\forall n$.111 To facilitate our analysis, we
ignore the non-linear coupling between the reflection amplitude and phase
shifts which has been discussed in [16, 3]. This practical model will be
considered in our future work.Because of the hardware limitation, IRS can only
choose its phase shifts from a finite number of discrete values. Therefore, we
assume the number of the quantization bits for the phase shifts at the IRS is
$B$, and the set of discrete phase shifts available at the IRS can be denoted
by
$\left\\{0,{\frac{2\pi}{2^{B}}},...,(2^{B}-1){\frac{2\pi}{2^{B}}}\right\\}$.
Then, the phase noise at IRS reflector $n$ can be written as
${\hat{\theta}_{n}}\sim{\mathcal{U}}\left[\frac{-\pi}{2^{B}},\frac{\pi}{2^{B}}\right],\forall
n$, where ${\mathcal{U}}$ denotes the uniform distribution. Therefore, the
phase shift matrix of the IRS with phase noise can be expressed as

$\displaystyle\tilde{\bf\Theta}={\rm{diag}}\left\\{{{\zeta_{1}}{e^{j{{\tilde{\theta}}_{1}}}},{\zeta_{2}}{e^{j{{\tilde{\theta}}_{2}}}},...,{\zeta_{N}}{e^{j{{\tilde{\theta}}_{N}}}}}\right\\},$
(2)

where ${\tilde{\theta}_{n}}={\theta_{n}}+{\hat{\theta}_{n}}$ denotes the
actual phase shift of the IRS element $n$ with $\theta_{n}$ being the designed
phase shift to be detailed later.

By properly deploying the IRS, the virtual LoS channel created by IRS in
mmWave systems can be strongly sparse. Thus, we adopt the widely used geometry
channel model [7, 8, 9] with only one dominant propagation path. The channel
from the user to the IRS and that from the IRS to the BS can be respectively
expressed as

$\displaystyle{\bf h_{1}}=\alpha{{\bf
a}_{N}}\left({\varphi_{r}^{a},\varphi_{r}^{e}}\right),{\bf H_{2}}=\beta{{\bf
a}_{M}}\left({\phi_{r}^{a},\phi_{r}^{e}}\right){\bf
a}_{N}^{H}\left({\varphi_{t}^{a},\varphi_{t}^{e}}\right),$ (3)

where $\alpha$ and $\beta$ are the link strengths. As shown in Fig. 1,
$\varphi_{r}^{a},\varphi_{r}^{e}$ ($\phi_{r}^{a},\phi_{r}^{e}$) are the
azimuth and elevation angles of arrival (AOA) at the IRS (BS),
$\varphi_{t}^{a},\varphi_{t}^{e}$ are the azimuth and elevation angles of
departure (AoD) at the IRS. ${{\bf
a}_{X}}\left({\vartheta^{a},\vartheta^{e}}\right)$ is the array response
vector of the USPA with size $\sqrt{X}\times\sqrt{X}$, which could be written
as

$\displaystyle\begin{array}[]{l}{{\bf
a}_{X}}\left({\vartheta^{a},\vartheta^{e}}\right)\\!=[1,...,{e^{j2\pi\frac{d}{\lambda}\left({x\sin\vartheta^{e}\sin\vartheta^{a}+y\sin\vartheta^{e}\cos\vartheta^{a}}\right)}},\\\
...,{e^{j2\pi\frac{d}{\lambda}\left({\left({\sqrt{X}-1}\right)\sin\vartheta^{e}\sin\vartheta^{a}+\left({\sqrt{X}-1}\right)\sin\vartheta^{e}\cos\vartheta^{a}}\right)}}{]^{T}},\end{array}$
(6)

where $\lambda$ and $d$ are the carrier wavelength and element spacing,
respectively. ${\rm{0}}\leq x,y\leq\sqrt{X}-1$ are the element indices of the
USPA. In this paper, we assume that angles in (3) can be perfectly estimated
based on some channel estimation methods[17, 18] and the case with imperfect
angular information as in [19] will be left for our future work.

To reduce the power consumption and hardware cost in mmWave band, low-
resolution ADCs are adopted at the BS. For analytical tractability, we adopt
the additive quantization noise model (AQNM) to characterize the impacts of
the low-resolution ADCs at the BSs. Thus, the quantizer outputs of the
received signal with phase noise can be written as [20]

$\displaystyle{{\bf y}_{q}}=\gamma{\bf y}+{{\bf
n}_{q}}=\gamma\left({\sqrt{P}\tilde{\bf g}{s}+{\bf n}}\right)+{{\bf
n}_{q}}=\gamma\sqrt{P}\tilde{\bf g}{s}+\gamma{\bf n}+{{\bf n}_{q}},$ (7)

where $\tilde{\bf g}={\bf H_{2}}\tilde{\bf\Theta}{\bf h_{1}}$ represents the
$M\times 1$ cascaded user-IRS-BS channel with phase noise, ${s}$ denotes the
transmit symbol, ${\bf n}$ is a normalized additive white Gaussian noise
vector with zero mean and unit variance. Therefore, $P$ represents the
normalized transmit SNR. Besides, $\gamma=1-\rho$, where $\rho$ is the inverse
of the signal-to-quantization-noise ratio and ${\bf n}_{q}$ denotes the
additive Gaussian quantization noise vector which is uncorrelated with $\bf
y$. The covariance matrix of ${\bf n}_{q}$ is given by

$\displaystyle{{\bf R}_{{{\bf n}_{q}}{{\bf
n}_{q}}}}=\gamma\left({1-\gamma}\right){\rm diag}\left({P\tilde{\bf
g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right).$ (8)

Considering the non-uniform scalar minimum mean-square-error (MMSE) quantizer
of the Gaussian random variables, the values of $\rho$ corresponding to the
quantization bits $b$ can be found from Table I[20].

TABLE I: $\rho$ versus Quantization Bits $b$ $b$ | 1 | 2 | 3 | 4 | 5 | $\geq 6$
---|---|---|---|---|---|---
$\rho$ | 0.3634 | 0.1175 | 0.03454 | 0.009497 | 0.002499 | $\frac{{\sqrt{3}\pi}}{2}{2^{-2b}}$

##  III Analysis of Achievable Uplink Rate

In this section, we analyze the uplink achievable rate and investigate the
effect of low-resolution ADCs at the BS and phase noise at the IRS. In
addition, we characterize the performance gap between our results and the
ideal case without quantization error or phase noise.

Consider the maximal-ratio-combining (MRC) receiver. Since the phase noise is
unknown, the receive beamforming is designed based on the available channel
state information (CSI) and on the optimized phase shift of the IRS. Thus, the
quantized signal received by the BS can be processed as

$\displaystyle{r}={{\bf g}^{H}}{{\bf y}_{q}}=\gamma\sqrt{P}{{\bf
g}^{H}}\tilde{\bf g}{s}+\gamma{{\bf g}^{H}}{\bf n}+{{\bf g}^{H}}{{\bf
n}_{q}},$ (9)

where ${\bf g}={\bf H_{2}}{\bf\Theta}{\bf h_{1}}$.

Then the uplink achievable rate can be written as

$\displaystyle R={\rm log}_{2}\\!\left(1+\frac{{{\gamma^{2}}P{{\left|{{{\bf
g}^{H}}\tilde{\bf g}}\right|}^{2}}}}{{{\gamma^{2}}{{\left\|{\bf
g}\right\|}^{2}}+\gamma\left({1-\gamma}\right){\bf g}^{H}{\rm
diag}\left({P\tilde{\bf g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right){\bf
g}}}\right).$ (10)

###  III-A Phase Shift Matrix Design

Since the phase noise at the IRS is unknown, we maximize the SNR by designing
the phase shifts without considering the existence of phase noise. Although
this design is sub-optimal, it has lower complexity in practical systems and
can lead to a tractable analytical expression. The optimization problem can be
formulated as follows

$\displaystyle\begin{array}[]{l}\max\limits_{\bf\Theta}\quad{\rm
SNR}^{ideal}=\frac{\gamma^{2}P\left|\mathbf{g}^{H}\mathbf{g}\right|^{2}}{\gamma^{2}\|\mathbf{g}\|^{2}+\gamma(1-\gamma)\mathbf{g}^{H}\operatorname{diag}\left(P\mathbf{gg}^{H}+\mathbf{I}\right)\mathbf{g}},\\\
\;{\rm s.t.}\;\;\quad{\theta}_{n}\in[0,2\pi),\forall n.\end{array}$ (13)

Recalling (3), ideal cascaded channel $\bf g$ can be rewritten as

$\displaystyle\begin{array}[]{l}\mathbf{g}=\mathbf{H}_{2}\mathbf{\Theta}\mathbf{h}_{1}=\alpha\beta\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)f({\bf\Theta}),\\\
\left|{\rm
g}_{i}\right|=\left|{\alpha\beta}\right|\left|f({\mathbf{\Theta}})\right|,\end{array}$
(16)

where
$f(\mathbf{\Theta})\triangleq\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\mathbf{\Theta}\mathbf{a}_{N}\left(\varphi_{r}^{a},\varphi_{r}^{e}\right)$
and $\mathrm{g}_{i}$ is the $i$-th entry of vector $\mathbf{g}$. Then, we can
respectively calculate the three terms in $\mathrm{SNR}^{ideal}$ as follows

$\displaystyle\begin{array}[]{l}\|\mathbf{g}\|^{2}=|\alpha\beta|^{2}\left\|\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\right\|^{2}|f(\mathbf{\Theta})|^{2}=|\alpha\beta|^{2}M|f(\mathbf{\Theta})|^{2},\\\
\left|\mathbf{g}^{H}\mathbf{g}\right|^{2}=\left(\|\mathbf{g}\|^{2}\right)^{2}=|\alpha\beta|^{4}M^{2}|f(\mathbf{\Theta})|^{4},\\\
\mathbf{g}^{H}\operatorname{diag}\left(P{\mathbf{g}}\mathbf{g}^{H}+\mathbf{I}\right)\mathbf{g}=\sum_{i=1}^{M}{\mathrm{g}^{*}_{i}}\left(P{\mathrm{g}_{i}}{\mathrm{g}^{*}_{i}}+1\right){\mathrm{g}_{i}}\\\
=\\!\\!\sum_{i=1}^{M}\\!\\!\left(\\!P\left|{\mathrm{g}_{i}}\right|^{4}\\!\\!+\\!\left|{\mathrm{g}_{i}}\right|^{2}\\!\right)\\!\\!=\\!M\\!\left(\\!P|\alpha\beta|^{4}|f(\mathbf{\Theta})|^{4}\\!+\\!|\alpha\beta|^{2}|f(\mathbf{\Theta})|^{2}\right).\end{array}$
(21)

By substituting (21) into $\mathrm{SNR}^{ideal}$, we can see that the maximal
$\mathrm{SNR}^{ideal}$ is achieved under maximal $|f(\mathbf{\Theta})|$. Thus,
we have

$\displaystyle\begin{array}[]{l}{{\bf\Theta}^{opt}}\\!=\\!\mathop{\arg\max}\limits_{\bf\Theta}|f(\mathbf{\Theta})|\\!=\\!\mathop{\arg\max}\limits_{\bf\Theta}{\left|{\sum_{n=1}^{N}{{e^{j2\pi\frac{d}{\lambda}\left({xp+yq}\right)+j{\theta_{n}}}}}}\\!\right|},\end{array}$
(23)

where $x=\lfloor{\left({n-1}\right)/\sqrt{N}}\rfloor$ and
$y=\left({n-1}\right)\bmod\sqrt{N}$. $\left\lfloor n\right\rfloor$ means
rounding $n$ toward the negative infinity and $\bmod$ means taking the
remainder after division.
$p=\sin\varphi_{r}^{e}\sin\varphi_{r}^{a}-\sin\varphi_{t}^{e}\sin\varphi_{t}^{a}$
and
$q=\sin\varphi_{r}^{e}\cos\varphi_{r}^{a}-\sin\varphi_{t}^{e}\cos\varphi_{t}^{a}$.
Based on the triangle inequality, the optimal phase shifts are given by

$\displaystyle\theta_{n}^{opt}=-2\pi\frac{d}{\lambda}\left({xp+yq}\right),\forall
n.$ (24)

###  III-B Uplink Rate Analysis

Based on the optimal phase shift design in (24), we next derive the
theoretical expression of the uplink achievable rate. The result is presented
in the following theorem.

###### Theorem 1

Considering the impacts of the quantization error at the BSs and phase noise
at the IRS, the uplink achievable rate for the IRS-assisted mmWave system with
large $N$ is

$\displaystyle R\mathop{\to}\limits^{a.s.}{\log_{2}}\left({1+\frac{{\gamma
P{{\left|{\alpha\beta}\right|}^{2}}{N^{2}}M{{\rm
sinc}^{2}}\left(\frac{\pi}{2^{B}}\right)}}{{\gamma+\left({1-\gamma}\right)\left({1+P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}{{{\mathop{\rm
sinc}\nolimits}}^{\rm{2}}}\left(\frac{\pi}{2^{B}}\right)}\right)}}}\right).$
(25)

Proof: See Appendix A. $\Box$

Theorem 1 characterize the impacts of the number of IRS reflecting elements
$N$, the number of antennas $M$ and transmit power $P$ on the uplink
achievable rate when considering the quantization error at the BS and the
phase noise at the IRS. The result in (25) shows that the uplink achievable
rate scaling law is ${\mathcal{O}}(\log_{2}(M))$, which is not affected by the
phase noise at the IRS. Besides, the following asymptotic results are
presented to help us gain a better understanding of Theorem 1.

###### Remark 1

For the ideal ADC without quantization error, i.e., $\gamma\to 1$, the uplink
achievable rate reduces to

$\displaystyle
R\mathop{\to}\limits^{a.s.}{\log_{2}}\left({1+P{{\left|{\alpha\beta}\right|}^{2}}{N^{2}}M{\rm
sinc^{2}}\left(\frac{\pi}{2^{B}}\right)}\right),$ (26)

which shows that the achievable rate can increase with $N$ infinitely.
Meanwhile, with a large $N$, (26) can be approximated as

$\displaystyle
R\to{\log_{2}}\left({P{{\left|{\alpha\beta}\right|}^{2}}{N^{2}}M}\right)+{\log_{2}}\left({{{\rm
sinc}^{2}}\left(\frac{\pi}{2^{B}}\right)}\right).$ (27)

This indicates that the phase noise at the IRS will cause a constant rate loss
${\log_{2}}\left({{{\rm sinc}^{2}}\left(\frac{\pi}{2^{B}}\right)}\right)$ when
the number of reflecting elements of the IRS becomes large.

###### Remark 2

With the quantization error and phase noise, when $P\\!\\!\to\\!\\!\infty$ or
$N\\!\\!\to\\!\\!\infty$, the uplink achievable rate reduces to

$\displaystyle R\to{\log_{2}}\left({1+\frac{{\gamma
M}}{{\left({1-\gamma}\right)}}}\right),$ (28)

which is independent of the phase noise at the IRS. Therefore, it is feasible
to use simple hardware with high phase noise on the IRS when deploying a large
number of reflecting elements or using high transmission power. In addition,
(28) shows that the achievable rate is mainly limited by the resolution of
ADCs and cannot increase infinitely when $P$ or $N$ becomes infinite.

Next, we give the following proposition to provide the design guidance for an
uplink IRS-aided mmWave system with an acceptable performance degradation.

###### Proposition 1

To design an uplink IRS-assisted mmWave system with rate degradation $\delta$,
the quantization precision of the IRS and the ADCs should respectively satisfy

$\displaystyle B\geq{\log_{2}}\pi-{\log_{2}}{\rm
arcsinc}\sqrt{\frac{{\hat{a}}}{{P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}\left({{2^{\delta}}\gamma
M-\hat{a}+\hat{a}\gamma}\right)}}},$ (29) $\displaystyle{\rho\leq
1-f(\hat{s})=1-}\frac{{\left({\left({1+M\hat{s}}\right){2^{-\delta}}-1}\right)\left({1+\hat{s}}\right)}}{{\left({\left({1+M\hat{s}}\right){2^{-\delta}}-1+M}\right)\hat{s}}},$
(30)

where ${{\hat{a}}}=1+\frac{{\gamma
P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}M}}{{1+\left({1-\gamma}\right)P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}}}-{2^{\delta}}$,
$\hat{s}=P{\left|{\alpha\beta}\right|^{2}}{N^{\rm{2}}}{{\mathop{\rm
sinc}\nolimits}^{\rm{2}}}\left({\frac{\pi}{{{2^{B}}}}}\right)$, the mapping
from $\rho$ to $b$ is given in Table I and
$\left(1-f\left(\hat{s}\right)\right)$ is a decreasing function with respect
to $N$, $P$, $B$ and $M$.

Proof: See Appendix B. $\Box$

On one hand, Proposition 1 reveals that the precision of ADCs at BSs should
increase with $N$, $P$, $B$ and $M$ to guarantee the acceptable rate loss. For
$N\to\infty$ or $P\to\infty$, we have $\rho\leq 0$ based on (30). For the
increase of $M$, the decreasing speed of $\rho$ decreases based on (45) which
shows that the first-order derivative is approaching zero when $M$ is very
large.

On the other hand, Proposition 1 shows that $B\geq 0$ when $N\to\infty$ or
$P\to\infty$, which means that the precision of IRS phase shifts is
insignificant under large $N$ or $P$. Besides, when $\delta=0$, we have
$B\geq\infty$ and $\rho\leq 0$. These results, which are consistent with the
Remark 1 and Remark 2, verify the accuracy of Proposition 1.

###### Remark 3

Assume that the transmission power $P$ is scaled with $M$ and $N$ according to
$P=\frac{{{E_{u}}}}{{M{N^{2}}}}$, where $E_{u}$ is fixed. When $M\to\infty$,
we can obtain

$\displaystyle
R\mathop{\to}\limits^{a.s.}{\log_{2}}\left({1+\gamma{E_{u}}{{\left|{\alpha\beta}\right|}^{2}}{{\rm
sinc}^{2}}\left(\frac{\pi}{2^{B}}\right)}\right).$ (31)

This indicates that with the assistance of the IRS, the transmit power can be
scaled down proportionally to $1/(MN^{2})$ to achieve the same rate of a
single-antenna non-IRS-aided system with transmit power ${\gamma{E_{u}}{{\rm
sinc}^{2}}\left(\frac{\pi}{2^{B}}\right)}$ and channel strength $\alpha\beta$.

##  IV Numerical Results

In this section, we validate our analytical results through simulations.
Unless otherwise stated, our parameters are set as $N=64$, $M=16$, $B=1$,
$b=2$, $\alpha^{2}=10^{-3}d_{1}^{-2.2}$, $\beta^{2}=10^{-3}d_{2}^{-2.2}$ [6],
$d_{1}=10$ m, $d_{2}=40$ m, $P=\frac{20\,{\rm dbm}}{-80\,{\rm dbm}}$.
Theoretical results are calculated using (25) while simulation results are
obtained by averaging over 100 independent calculations of (10).

Figure 2: Uplink achievable rate versus $N$.

Figure 3: The impact of IRS quantization bits $B$.

Figure 4: Rate versus the number of antennas $M$.

Fig. 4 shows the the impacts of the number of IRS reflecting elements on the
uplink achievable rate under four different scenarios. It can be observed that
our analytical results match well with the simulation results, which verifies
the correctness of our derivations. Meanwhile, the simulation results of Fig.
4 are consistent with our discussions in Remark 1 and Remark 2. With an ideal
ADC, we can see that the rate can be increased infinitely by enlarging the
size of IRS, and the gap between the ideal IRS (i.e., without phase noise) and
non-ideal IRS tends to be a constant when $N$ becomes large. On the contrary,
when considering the effect of the low-resolution ADC, the rate with/without
phase noise finally approaches the same upper bound.

Fig. 4 shows the influence of the phase noise under different ADCs
resolutions. We can see that phase noise at the IRS has a larger impact when
the ADC has a higher resolution. When the precision of the ADC decreases, only
the serious phase noise, i.e., $B=1$, will bring obvious degradation of system
performance. However, even with ideal ADCs, a 3-bits quantized IRS is
sufficient to achieve the near-optimal rate performance. This finding
validates the feasibility of using low-quality hardware at IRS (e.g., $B=1,2$)
when the BSs are equipped with low-resolution ADCs.

Fig. 4 shows the uplink achievable rate scaling law with respect to the number
of the antennas at the BS. We can see that the achievable rate increases
logarithmically with $M$ when $M$ is large. Besides, we can also see that the
impact of phase noise becomes marginal when $N$ is large.

##  V Conclusion

In this paper, we have analyzed the rate performance of an uplink IRS-aided
MISO system with hardware limitations at both the BS and the IRS. Using our
derived expressions, we have shown that the rate scaling law is
${\mathcal{O}}(\log_{2}(M))$ which is not influenced by the phase noise at the
IRS. Our results have also shown that when the number of reflecting elements
of the IRS approaches infinity, the rate is only limited by the precision of
the ADCs. Numerical simulations have verified the accuracy of our analytical
results.

## Appendix A

To derive the expression of the rate in (10), we need to calculate ${{{\bf
g}^{H}}\tilde{\bf g}}$, ${{{\left\|{\bf g}\right\|}^{2}}}$ and ${\bf
g}^{H}{\rm diag}\left({P\tilde{\bf g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right){\bf
g}$, respectively. Please note that in these terms, the only random variable
is phase noise ${\hat{\theta}}_{n}$, which will converge to its expectation
when $N\to\infty$ due to the law of large number.

First, by using the optimal phase shift design in (24), the cascaded channels
with and without phase noise can be respectively expressed as

$\displaystyle\begin{array}[]{l}{\bf g}={\bf H_{2}}{\bf\Theta}{\bf
h_{1}}=\alpha\beta N{{\bf a}_{M}}\left({\phi_{r}^{a},\phi_{r}^{e}}\right),\\\
\tilde{\bf g}={\bf H_{2}}\tilde{\bf\Theta}{\bf h_{1}}=\alpha\beta{{\bf
a}_{M}}\left({\phi_{r}^{a},\phi_{r}^{e}}\right)\sum_{n=1}^{N}{{e^{j{{\hat{\theta}}_{n}}}}}.\end{array}$
(34)

Then using the result in (34), we have

$\displaystyle{\left\|{\bf
g}\right\|^{2}}=\mathbf{g}^{H}{\mathbf{g}}\\!=\\!|\alpha\beta|^{2}N^{2}{\left\|\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\right\|^{2}}\\!=\\!{\left|{\alpha\beta}\right|^{2}}{N^{2}}M,$
(35)
$\displaystyle\mathbf{g}^{H}\tilde{\mathbf{g}}=\\!|\alpha\beta|^{2}N{\left\|\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\right\|^{2}}\\!\sum\limits_{n=1}^{N}\\!e^{j\hat{\theta}_{n}}\\!=\\!|\alpha\beta|^{2}NM\\!\sum\limits_{n=1}^{N}\\!e^{j\hat{\theta}_{n}}.$
(36)

For large $N$, by using strong law of large numbers and the continuous mapping
theorem, we can obtain

$\displaystyle\frac{1}{N}\\!\sum_{n=1}^{N}\\!{{e^{\pm
j{{\hat{\theta}}_{n}}}}}\mathop{\to}\limits^{a.s.}\mathbb{E}\\!\left[{{e^{\pm
j{{\hat{\theta}}_{n}}}}}\right]\\!{\mathop{=}\limits^{\left(a\right)}}\mathbb{E}\\!\left[{\cos\\!\left({{{\hat{\theta}}_{n}}}\right)}\right]\\!{\mathop{=}\limits^{\left(b\right)}}{\mathop{\rm
sinc}\nolimits}\\!\left(\frac{\pi}{2^{B}}\\!\right),$ (37)

where $(a)$ applies the symmetry of the odd function
$\sin\left({{{\hat{\theta}}_{n}}}\right)$ with
${\hat{\theta}_{n}}\sim{\mathcal{U}}\left[\frac{-\pi}{2^{B}},\frac{\pi}{2^{B}}\right]$,
and $(b)$ is calculated based on the probability density function of
${\hat{\theta}}_{n}$. Therefore, we have

$\displaystyle{{{\bf g}^{H}}\tilde{\bf
g}}{\mathop{\to}\limits^{a.s.}}{{{\left|{\alpha\beta}\right|}^{2}}{N^{2}}{M}}{{{{\mathop{\rm
sinc}\nolimits}}}\left(\frac{\pi}{2^{B}}\right)}.$ (38)

Then we calculate the term ${\bf g}^{H}{\rm diag}\left({P\tilde{\bf
g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right){\bf g}$. The $m$th diagonal element of
${\rm diag}\left({P\tilde{\bf g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right)$ can be
written as

$\displaystyle{\left[{{\rm diag}\left({P\tilde{\bf g}{{\tilde{\bf
g}}^{H}}\\!+\\!{\bf
I}}\right)}\right]_{mm}}\\!=\\!1\\!+\\!P{\left|{\alpha\beta}\right|^{2}}\sum\limits_{n=1}^{N}{{e^{j{{\hat{\theta}}_{n}}}}}\\!\\!\sum\limits_{n=1}^{N}{{e^{-j{{\hat{\theta}}_{n}}}}}.$
(39)

Thus, by applying (37), we have

$\displaystyle{\bf g}^{H}{\rm diag}\left({P\tilde{\bf g}{{\tilde{\bf
g}}^{H}}+{\bf I}}\right){\bf
g}={\left|{\alpha\beta}\right|^{2}}{N^{2}}M\left[{{\rm diag}\left({P\tilde{\bf
g}{{\tilde{\bf g}}^{H}}+{\bf I}}\right)}\right]_{mm}$
$\displaystyle\mathop{\to}\limits^{a.s.}{\left|{\alpha\beta}\right|^{2}}{N^{2}}M\left({1+P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}{{{\mathop{\rm
sinc}\nolimits}}^{\rm{2}}}\left(\frac{\pi}{2^{B}}\right)}\right).$ (40)

Finally, substituting (35), (38) and (A) into (10) and with some simple
manipulations, we can complete the proof.

## Appendix B

In order to ensure that the performance loss is lower than $\delta$, for the
quantization bits of IRS (i.e., $B$) and the resolution of ADCs (i.e.,
$\gamma$), it should respectively satisfy the following inequalities

$\displaystyle\begin{array}[]{l}{\log_{2}}\left({1+\frac{{\gamma
P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}M}}{{\gamma+\left({1-\gamma}\right)\left({1+P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}}\right)}}}\right)-R\left({\gamma,B}\right)\leq\delta,\\\
{\log_{2}}\left({1+P{{\left|{\alpha\beta}\right|}^{2}}{N^{\rm{2}}}M{{{\mathop{\rm
sinc}\nolimits}}^{\rm{2}}}\\!\left({\frac{\pi}{{{2^{B}}}}}\right)}\right)-R\left({\gamma,B}\right)\leq\delta,\end{array}$
(43)

where $R\left({\gamma,B}\right)$ is the general rate expression (25) presented
in Theorem 1.

By solving inequalities in (43), we can obtain the constraints of $B$ and
$\rho$ in (29) and (30), respectively.

Next, based on the expression of $f(\hat{s})$, the first-order derivative of
$f(\hat{s})$ with respect to $\hat{s}$ and $M$ can be respectively written as

$\displaystyle\frac{{\partial
f\left({\hat{s}}\right)}}{{\partial\hat{s}}}\mathop{=}\limits^{\left(c\right)}\frac{{2M\hat{s}{2^{-\delta}}\left({1-{2^{-\delta}}}\right)+\left({M-1}\right)\left({1-{2^{-\delta}}}\right)}}{{{{\left({\hat{s}{2^{-\delta}}+M{{\hat{s}}^{2}}{2^{-\delta}}-\hat{s}+M\hat{s}}\right)}^{2}}}}$
$\displaystyle\qquad\quad+\frac{{{M^{2}}{{\hat{s}}^{2}}{2^{-\delta}}\left({1-{2^{-\delta}}}\right)+{2^{-\delta}}\left({1-{2^{-\delta}}}\right)}}{{{{\left({\hat{s}{2^{-\delta}}+M{{\hat{s}}^{2}}{2^{-\delta}}-\hat{s}+M\hat{s}}\right)}^{2}}}}\geq
0,$ (44) $\displaystyle\frac{{\partial f\left({\hat{s}}\right)}}{{\partial
M}}\mathop{=}\limits^{\left(c\right)}\frac{{\hat{s}\left({1-{2^{-\delta}}}\right)+{{\hat{s}}^{2}}\left({1-{2^{-\delta}}}\right)}}{{{{\left({\hat{s}{2^{-\delta}}+M{{\hat{s}}^{2}}{2^{-\delta}}-\hat{s}+M\hat{s}}\right)}^{2}}}}\geq
0,$ (45)

where $(c)$ is calculated based on the derivative law of fractions and through
some algebraic simplifications.

Since $\hat{s}$ is increasing with $N$, $P$ and $B$, thus we can complete the
proof.

## References

  * [1] M. Xiao _et al._ , “Millimeter wave communications for future mobile networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 9, pp. 1909–1935, Sep. 2017.
  * [2] M. Di Renzo _et al._ , “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2450–2525, Nov. 2020.
  * [3] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless communications: A tutorial,” _arXiv preprint arXiv:2007.02759_ , 2020.
  * [4] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5218–5233, Aug 2020.
  * [5] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [6] C. Pan _et al._ , “Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 8, pp. 1719–1734, Aug. 2020.
  * [7] G. Zhou, C. Pan, H. Ren, K. Wang, M. Elkashlan, and M. Di Renzo, “Stochastic learning-based robust beamforming design for ris-aided millimeter-wave systems in the presence of random blockages,” _arXiv preprint arXiv:2009.09716_ , 2020.
  * [8] P. Wang, J. Fang, X. Yuan, Z. Chen, and H. Li, “Intelligent reflecting surface-assisted millimeter wave communications: Joint active and passive precoding design,” _IEEE Transactions on Vehicular Technology_ , eatly access 2020.
  * [9] C. Pradhan, A. Li, L. Song, B. Vucetic, and Y. Li, “Hybrid precoding design for reconfigurable intelligent surface aided mmwave communication systems,” _IEEE Wireless Communications Letters_ , vol. 9, no. 7, pp. 1041–1045, Jul. 2020.
  * [10] A. Alkhateeb, J. Mo, N. Gonzalez-Prelcic, and R. W. Heath, “MIMO precoding and combining solutions for millimeter-wave systems,” _IEEE Communications Magazine_ , vol. 52, no. 12, pp. 122–131, Dec. 2014.
  * [11] M. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” _IEEE Wireless Communications Letters_ , vol. 9, no. 2, pp. 184–188, Feb. 2019.
  * [12] S. Zhou, W. Xu, K. Wang, M. Di Renzo, and M.-S. Alouini, “Spectral and energy efficiency of IRS-assisted MISO communication with hardware impairments,” _IEEE Wireless Communications Letters_ , 2020.
  * [13] Y. Han, W. Tang, S. Jin, C.-K. Wen, and X. Ma, “Large intelligent surface-assisted wireless communication exploiting statistical CSI,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 8, pp. 8238–8242, Aug. 2019.
  * [14] D. Li, “Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors,” _IEEE Communications Letters_ , vol. 24, no. 8, pp. 1646–1650, Aug. 2020.
  * [15] Z. Xing and R. Wang, “Achievable rate analyses and phase shift optimizations on intelligent reflecting surface with hardware impairments,” _arXiv preprint arXiv:2005.14411_ , 2020.
  * [16] S. Abeywickrama, R. Zhang, Q. Wu, and C. Yuen, “Intelligent reflecting surface: Practical phase shift model and beamforming optimization,” _IEEE Transactions on Communications_ , vol. 68, no. 9, pp. 5849–5863, Sep. 2020.
  * [17] W. Wang and W. Zhang, “Joint beam training and positioning for intelligent reflecting surfaces assisted millimeter wave communications,” _arXiv preprint arXiv:2009.03536_ , 2020.
  * [18] B. Deepak, R. Sankar, and S. P. Chepuri, “Channel estimation in reconfigurable intelligent surface assisted mmWave MIMO systems,” _arXiv preprint arXiv:2011.00900_ , 2020.
  * [19] X. Hu, C. Zhong, M. Alouini, and Z. Zhang, “Robust design for IRS-aided communication systems with user location uncertainty,” _IEEE Wireless Communications Letters_ , early access 2020.
  * [20] L. Fan, S. Jin, C.-K. Wen, and H. Zhang, “Uplink achievable rate for massive MIMO systems with low-resolution ADC,” _IEEE Communications Letters_ , vol. 19, no. 12, pp. 2186–2189, Dec. 2015.


# The Sensitivity of Sea Ice Brine Fraction to the Freezing Temperature and
Orientation

Kial D. STEWART  1,2 William PALM  1,2 Callum J. SHAKESPEARE  1,2 Noa
KRAITZMAN3 1Climate & Fluid Physics Laboratory, Australian National University
2Australian Research Council Centre of Excellence for Climate Extremes
3Macquarie University
Correspondence: Kial Stewart<EMAIL_ADDRESS>
###### Abstract

Pound for pound, sea ice is the most important component of Earth’s climate
system. The changing conditions in which sea ice forms and exists are likely
to affect the properties of sea ice itself, and potential climate feedbacks
need to be identified and understood to improve future projections. Here we
perform a set of idealised experiments which model sea ice growth under a
range of ambient salinities, freezing temperatures, and freezing orientations.
The results confirm existing theories; sea ice growth rate is largest for
cooler freezing temperatures, fresher ambient salinities, and bottom oriented
freezing configuration. Our primary metric of interest is the brine mass
fraction (the mass ratio of brine inclusions to the total sea ice), which we
quantify and determine the sensitivity of with respect to changes in freezing
conditions. We find that the brine mass fraction of our model sea ice is most
sensitive to freezing temperature, and increases 2.5% per 1∘C increase of
freezing temperature. This finding suggests that future sea ice in warmer
climate states will retain more brine, with subsequent flow on effects for
circulations driven by brine rejection in the high-latitude oceans.

##  1 Introduction

Sea ice is a vital component of Earth’s climate system, covering more than 10%
of the global oceans at least part of each year (??). Its high albedo reflects
between 80-90% of the incoming solar radiation, thereby strongly limiting the
radiative heating of polar surface waters. Brine rejection from sea ice
increases the salinity of high-latitude surface ocean waters that ultimately
sink and play a key role in the global thermohaline overturning circulation.
Sea ice is also highly sensitive to changes in climate, and is currently
exhibiting repeated record low areal coverage (e.g., ?), which serves as an
important measure of how the Earth system is responding to climate change.

Despite the importance of its role in Earth’s climate, and its unique set of
dynamics and thermodynamics, the representation of sea ice in coupled climate
models tends to mirror that of the oceanic circulation models. That is, the
sea ice field is discretised onto the ocean’s horizontal grid, the sea ice
state and characteristics are represented with bulk values, and the
dynamics/thermodynamics are approximated with a set of continuous partial
differential equations that are solved in time. Theories and limited existing
observations are used to develop parameterisations of small scale dynamics
that relate sub-grid scale processes to the larger scales that are resolved by
the horizontal grid. Many dynamics specific to sea ice, in particular its
rheology, need to be omitted because they are not well represented by the
ocean’s horizontal grid. Nevertheless, sea ice model development continues to
improve and include additional processes, whose representation in turn often
require new parameterisations. For example, the latest version of the CICE sea
ice model (CICE6 v6.1.2) includes a capability for a floe size distribution,
which then permits lateral sea ice growth/melt dynamics, which requires
specific parameterisations that are distinct from vertical sea ice
growth/melt. This new feature raises the question of whether the direction of
sea ice growth changes the bulk characteristics of the sea ice itself (e.g.,
?).

The difficulties associated with the realistic representation of sea ice in
climate models are a major scientific challenge for predicting future climate
states. Sea ice is a particularly sensitive entity in the climate system and
prone to feedback cycles, as demonstrated by the ice-albedo phenomenon (e.g.,
?). The polar amplification of global warming (e.g., ?) will tend to the warm
surface atmosphere at higher latitudes faster than lower latitudes; this will
result in sea ice growing in conditions that are relatively warmer than
before. Understanding how sea ice properties respond to warmer freezing
temperatures is an important step for improving sea ice and climate models.

Here, we investigate the impacts of freezing temperature and orientation on
the brine mass fraction of sea ice by conducting a suite of laboratory
experiments to model sea ice growth across a range of freezing temperatures
and orientations. The bulk brine mass fraction of the model sea ice is
estimated using two independent approaches, and found to be sensitive to both
the freezing temperature and orientation across a wide range of ambient
salinities. In §2 we provide background context for our study and subsequent
analysis. In §3 we describe the laboratory apparatus, methodology and
analysis, and present and discuss the results of the experiments in §4. We
provide a brief commentary on the geophysical implications of our results in
§5.

##  2 Background

Seawater is a multi-component liquid wherein its constituents have different
freezing points. In such a system, the freezing process typically results in
the formation of a porous mixture of solid and liquid inclusions known as a
mushy layer (?). Solidification at temperatures above the melting temperature
of the constituents results in the solvent freezing whilst the solutes remain
liquid. In the case of an aqueous solution containing salts, like seawater, as
solidification occurs, finger-like solid crystals known as dendrites grow in
advance of the solidifying face. These dendrites cannot incorporate the
dissolved salts into their matrix and expel it into the solution between the
solid crystals (?). This enriched saline solution is referred to as brine. As
the dendrites continue to grow and the ice face advances, some of the liquid
brine becomes trapped within the solid ice, forming brine inclusions. The
mixture of solid ice and brine inclusions forms the mushy layer (?). This
mixture is generally classified by its solid mass fraction, which is the ratio
of the solid ice mass to the total mass of the mush; the solid mass fraction
is related to the brine mass fraction (the ratio of the liquid brine mass to
the total mushy layer mass, which is often employed in sea ice models) as they
sum to unity.

Mushy layers are prevalent on Earth and form from both natural and industrial
processes, such as the solidification of Earth’s core, magma, and metallic
alloys (???). In the case of seawater, during the formation of sea ice, salts
from the seawater are rejected from the ice lattice to form highly saline
brine inclusions. The high salinity of these brine inclusions results in them
remaining in liquid phase within the ice; that is, in temperature–salinity
space they exist to the right of the of the liquidus curve (e.g., Fig. 1).

Mushy layers are subject to chemical and thermodynamic processes, and as such
they react and evolve as the mushy layer grows (??). For example, brine
inclusions must be in thermodynamic equilibrium within the ice. If they are
not, heat will flow to negate the temperature gradient such that thermodynamic
equilibrium is established. Therefore, due to the low temperatures within
mushy layers, in order for the brine inclusions to remain liquid, they must
exist on their liquidus curve, the solid-liquid phase transition point of a
given substance. This requires brine inclusions to vary their chemical and
temperature characteristics as the mushy layer varies in temperature, which
results in a reactive medium.

? demonstrated that solidification of solutions with multiple dissolved salts
below the eutectic temperature (the minimum temperature where all salts are
able to be in solution) creates a compositional mushy layer; that is, an ice
matrix which contains crystals of solid salt. Here we avoid the added
complexity of compositional mushy layers and focus on solidification at
temperatures above the eutectic temperature for our aqueous solution. We model
the seawater and sea ice with a single-component liquid that is sodium
chloride (NaCl) and fresh tapwater, which has a single eutectic temperature of
$T_{e}=-21.1^{\circ}$C. Note that actual seawater and sea ice, as a material
composed of several different salts, does not have a eutectic point, and there
is still liquid brine for temperatures below $-70^{\circ}$C, however
precipitation begins at $T_{e}<-21.1^{\circ}$C (e.g., ?), and modern sea ice
models employ an equation of state with an imposed a eutectic temperature of
$-36.2^{\circ}$C (e.g., ?).

The temperature difference between the relatively warmer ocean and relatively
cooler atmosphere means that vertical temperature gradients exist within sea
ice. It follows that a similar temperature gradient exists within and
surrounding individual brine inclusions (e.g., ?). For a brine inclusion to
maintain phase equilibrium, a salinity gradient is necessary within the liquid
brine inclusions, which requires saltier, cooler liquid to exist at the top of
the inclusion (?). Obviously this arrangement is gravitationally unstable, and
the brine inclusion itself is likely to be well-mixed. That is, while the
bulk-average of the brine inclusion is likely to be in thermal and chemical
equilibrium, the upper and lower regions of the brine inclusion are not able
to achieve equilibrium with their respective surrounds; the sea ice
surrounding the upper regions of the inclusion will tend to freeze, while sea
ice surrounding the lower regions will tend to melt (e.g., ?). This process
allows the brine inclusions to tunnel through the sea ice which results in the
brine inclusions travelling from the upper, cooler to the lower, warmer
regions of the sea ice, and ultimately into the ocean below. Without a source
of salt in the upper sea ice region, the tunnelling of brine results in a
monotonic decrease in the bulk sea ice salinity in time. Additionally, the
direction of the temperature gradient within the sea ice determines the
tunnelling direction of brine inclusions; that is, horizontal gradients should
result in horizontal tunnels.

The amalgamation of tunnelling brine inclusions can develop solid free
chimneys referred to as brine channels (?). Brine channels provide a link
between the internals of the ice and the surrounding solution. In the case of
sea ice, the highly saline and near-freezing brine inclusions are much denser
than the ocean waters immediately below. Therefore, brine inclusions are
negatively-buoyant and brine channels can facilitate convective fluid
processes between sea ice and sea water (?). Convective motion expels saline
brine into the underlying ocean, in a process known as gravity drainage. As
brine channels provide a direct and active mechanism to drive flow between sea
ice and the ocean, gravity drainage governs the rate of brine rejection. Due
to continuity, the volume of brine expelled from ice must be replaced,
resulting in relatively less saline seawater replacing brine in the mushy
layer (?). As thermodynamic and chemical equilibrium must again be
established, the less saline seawater in the ice freezes, re-initiating
solidification and hence, increasing the solid mass fraction, and altering the
mushy layer structure. Increasing the solid mass fraction of ice decreases its
permeability and thus a point is reached where convective motions can no
longer be sustained, halting gravity drainage. Therefore, gravity drainage is
a mechanism for brine rejection that is most active during the early stages of
sea ice growth. Orientation of the freezing sea ice face will likely determine
how effective gravity drainage is in facilitating brine rejection.

Figure 1: A partial phase diagram for aqueous sodium chloride solutions. The
blue region is cooler than the eutectic temperature $T_{e}=-21.1^{\circ}$C, at
which sodium chloride crystal precipitate out of solution. The white region
represents the mushy layer regime, and the yellow region represents the regime
of aqueous solution wherein the salinity is too large for ice to form for the
given temperature; these regimes are separated by the liquidus curve, which is
given by $T_{liq}$ (Eqn. 8) for temperatures warmer than eutectic temperature.
The crosses, circles, and triangles represent the initial ambient salinities
and freezing plate temperatures of the experiments with top, side, and bottom
freezing plate orientations, respectively.

##  3 Laboratory Experiments

###  3.1 Apparatus

Figure 2: A two-dimensional schematic of the experiment apparatus.

Our experiments investigating ice growth were conducted in the Climate & Fluid
Physics (CFP) laboratory at the Australian National University. The
experimental apparatus consists of a well-insulated, sealed, fluid-filled tank
with internal dimensions $30.5{\times}21.5{\times}17.5$ cm3 (see Figure 2), of
a similar design to that used by ?. The two smallest area sidewalls of the
tank are copper plates that are in direct thermal contact with controllable
heat exchangers that allow the sidewall temperatures to be independently
prescribed; the other four sidewalls are double-glazed perspex. The tank is
positioned between an illuminated white LED screen and a Basler AG camera that
views the entire tank volume and configured such that it is able to clearly
distinguish between the ambient fluid and the growing ice. The apparatus is in
a section of the CFP laboratory that is climate controlled and maintained at a
constant temperature of $22.5\pm 0.2^{\circ}C$.

Two independent Julabo FP50 HL refrigerating/heating constant temperature
baths are used to prescribe the copper sidewall temperatures. These constant
temperature Julabo units pump ethanol through the sidewall heat exchangers,
thereby providing a operational sidewall temperature range that easily spans
from below the eutectic temperature of salty water ($T_{e}=-21.1^{\circ}$C) to
above the freezing point of fresh water. During operation, one of these copper
sidewall plates is set to a temperature below the freezing point of the fluid,
with the other plate used to maintain the temperature of the ambient fluid;
these are referred to as the “freezing” and “ambient” plates, respectively. So
while the freezing plate initiates ice growth, the primary purpose of the
ambient plate is to negate any unwanted heat flux into the tank from the
laboratory, thereby maintaining the temperature of the ambient close to its
freezing temperature. The ambient fluid is kept well-mixed and homogeneous by
an external pump that gently recirculates the ambient in the vicinity of the
ambient plate and inhibits the formation of stratification. The temperature of
the ambient is logged using a thermistor positioned in the centre of the tank.

The effect of the freezing direction is one of the variables of interest. For
this, the tank is mounted in a semicircular frame that is able to be tilted
such that the freezing plate can be oriented at the top, side, or bottom of
the tank.

As ice grows against the freezing plate, there is an increase in total volume
that is associated with the phase change from liquid to ice. This volume
increase is accommodated by a displacement of ambient fluid through an
expansion tube and into a beaker that is on a scale, which records the
measured mass in time.

###  3.2 Methodology

The tank was tilted in the desired orientation and the tank, expansion tube,
and scale beaker were carefully filled with a de-aired solution of water with
a precisely known salinity. The open end of the expansion tube was placed
within the scale beaker below the water level, thereby creating a closed
system that is able to accommodate any volume changes in the tank. The mixing
pump was switched on to keep the ambient fluid well-mixed. Both of the
constant temperature sidewalls were set to approximately 2∘C above the
estimated freezing temperature for the particular ambient salinity, and the
system was left to come to thermal equilibrium, which typically involved a
slight contraction of volume in the tank and an associated drawing of fluid
from the scale beaker. Over time, the mass measured by the scale adjusts to
become steady, at which time the experiment is ready to begin.

The freezing plate temperature is lowered to a desired point below the ambient
freezing temperature, initiating ice growth on the plate. The ambient plate
temperature is maintained just above the ambient freezing temperature. These
thermal boundary conditions were maintained for a period of over 24 hours,
during which time the ambient temperature and beaker mass were logged, and
high-resolution images from the Basler AG camera recorded at 30 minute
intervals. After this $\sim$24 hour logging period, the freezing plate
temperature was returned to that of the ambient plate and the ice was allowed
to melt back into the mixed ambient fluid. Once the ambient fluid returns to
thermal equilibrium, the next experiment is commenced with a different
prescribed freezing plate temperature and/or tank orientation. The ambient
fluid was completely replaced between experiments with different initial
ambient salinities.

###  3.3 Parameter Space

The freezing plate temperatures $T_{FP}$ explored here are $-10^{\circ}$C,
$-15^{\circ}$C, and $-20^{\circ}$C, which are all warmer than the eutectic
temperature of salty water at which compositional mushy layers form
($T_{e}=-21.1^{\circ}$C). The initial ambient salinities employed here are
approximately 17, 33, 48, 73, and 102 g/kg. For initial salinities greater
than 45 g/kg, the warmest freezing temperature setting ($T_{FP}=-10^{\circ}$C)
does not produce sufficiently detectable ice, so only freezing temperatures of
$T_{FP}=-15^{\circ}$C and $T_{FP}=-20^{\circ}$C are used. For the highest
salinity experiments (102 g/kg), ice produced in the bottom orientation was
too buoyant and would detach from the freezing plate and float up through the
ambient and melt, so only the top and side orientations are included here.
This approach provides a set of 34 different experiments, for which the
initial ambient salinities, freezing temperature settings, and tank
orientations are depicted in Figure 1.

###  3.4 Analysis: Calculating the Bulk Solid Mass Fraction

The quantitative diagnostics here are the logs of the ambient fluid
temperature and beaker mass, and the high-resolution images, which combine to
provide estimates of the bulk solid mass fraction $\phi$ of the mushy layer.
The solid mass fraction is defined by the ratio of the mass of solid ice to
the total mass of the mushy layer,

$\phi=\frac{m_{s}}{m_{b}+m_{s}},$ (1)

where $m_{s}$ and $m_{b}$ refer to mass of solid ice and liquid brine (kg),
respectively. A solid mass fraction of $\phi=1$ represents solid ice with no
liquid brine, and a value of $\phi=0$ is entirely liquid brine without any
solid ice. Note that the solid mass fraction $\phi$ and the brine mass
fraction $\phi_{l}$ ($\phi_{l}=m_{b}/\left(m_{b}+m_{s}\right)$), which is
often used in sea ice models in place of the solid mass fraction (e.g., ?),
sum to unity; $\phi+\phi_{l}=1$. The definition above describes a bulk solid
mass fraction which is useful as it provides a domain average of the varying
fine scale structures, such that the large-scale properties of the mushy layer
can be characterised. Here, we follow ? by employing two distinct approaches
to estimate the bulk solid mass fraction in our experiments; these are based
on the conservation of salt and the conservation of mass, respectively, and
are described below. Note that in reality the solid mass fraction exhibits
variability within the sea ice, which is not able to be captured with our bulk
conservation methods employed here.

####  3.4.1 Bulk Salt Conservation Method

The conservation of salt requires the total amount of salt in the system to
remain constant throughout the experiment, that is,

$\underbrace{S_{o}\rho_{o}V_{t}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Initial
Salt}}=\underbrace{(1-\phi)S_{b}\rho_{b}V_{m}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Brine
Inclusions}}+\underbrace{S_{a}\rho_{a}(V_{t}-V_{m})\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Ambient
Fluid}}+\underbrace{S_{e}m_{e}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Expelled
Fluid}},$ (2)

where $S_{o}$, $S_{b}$, $S_{a}$, and $S_{e}$ refer to the salinities of the
initial ambient fluid, brine inclusions, ambient fluid during ice growth, and
expelled fluid in the beaker, respectively (g/kg); $\rho_{o}$, $\rho_{b}$, and
$\rho_{a}$ are the densities of the initial ambient fluid, brine inclusions,
and ambient fluid during ice growth, respectively (kg/m3); $V_{t}$ and $V_{m}$
refer to the volumes of the tank and mushy layer, respectively (m3); and
$m_{e}$ is the mass of ambient fluid expelled into the beaker (kg).

Rearranging Equation 2 for the solid mass fraction produces,

$\phi_{S}=\frac{S_{o}\rho_{o}V_{t}-S_{b}\rho_{b}V_{m}-S_{a}\rho_{a}(V_{t}-V_{m})-S_{e}m_{e}}{-S_{b}\rho_{b}V_{m}},$
(3)

where we use the subscript $S$ identify this solid mass fraction estimate as
that provided by the salt conservation method.

####  3.4.2 Bulk Mass Conservation Method

The conservation of mass provides a mass balance model for the solid mass
fraction $\phi$ where the mass of the expelled fluid is directly measured, and
the masses of the initial system, solid ice, brine inclusions, evolving
ambient fluid are estimated by the products of their respective densities and
volumes. That is,

$\underbrace{\rho_{o}V_{t}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Initial
Mass}}=\underbrace{\phi\rho_{ice}V_{m}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Solid
Ice}}+\underbrace{(1-\phi)\rho_{b}V_{m}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Brine
Inclusions}}+\underbrace{\rho_{a}(V_{t}-V_{m})\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Ambient
Fluid}}+\underbrace{m_{e}\rule[-2.0pt]{0.0pt}{5.0pt}}_{\mbox{Expelled
Fluid}},$ (4)

where $\rho_{o}$, $\rho_{ice}$, $\rho_{a}$ and $\rho_{b}$ are the densities of
the initial ambient fluid, solid ice, ambient fluid during ice growth, and
brine channel fluid, respectively (kg/m3); $V_{t}$ and $V_{m}$ refer to the
volumes of the tank and mushy layer, respectively (m3); and $m_{e}$ is the
mass of ambient fluid expelled into the beaker (kg). Equation 4 can then be
rearranged for $\phi$, to provide the mass conservation method for estimating
the solid mass fraction of mushy layer, that is,

$\phi_{M}=\frac{\rho_{o}V_{t}-m_{e}-\rho_{a}(V_{t}-V_{m})-\rho_{b}V_{m}}{(\rho_{ice}-\rho_{b})V_{m}},$
(5)

where we use the subscript $M$ to identify this solid mass fraction estimate
as that provided by the mass conservation method.

####  3.4.3 Quantifying Terms for the Conservation Methods

Table 1: Terms needed for the Conservation Methods Term | Units | Description | Provenance
---|---|---|---
$h$ | m | mushy layer thickness | estimated from photos
$h_{o}$ | m | initial mushy layer thickness | estimated from photos
$h_{f}$ | m | final mushy layer thickness | estimated from photos
$V_{t}$ | m3 | volume of tank | directly measured
$V_{m}$ | m3 | volume of mushy layer | needs $h$
$m_{e}$ | kg | mass of expelled fluid | directly measured and logged
$\rho_{e}$ | kg/m3 | density of expelled fluid | directly measured
$S_{e}$ | g/kg | salinity of expelled fluid | Eqn. 6; needs $\rho_{e}$
$\rho_{o}$ | kg/m3 | initial ambient fluid density | directly measured
$S_{o}$ | g/kg | initial ambient fluid salinity | Eqn. 6; needs $\rho_{o}$
$\rho_{ice}$ | kg/m3 | density of solid ice | Eqn. 7; needs $T_{ice}$
$T_{ice}$ |  ∘C | interior temperature of solid ice | Eqn. 9; needs $T_{liq}$
$T_{liq}$ |  ∘C | liquidus temperature | Eqn. 8; needs $S_{a}$
$S_{b}$ | g/kg | salinity of brine inclusions | Eqn. 10; needs $T_{ice}$
$\rho_{b}$ | kg/m3 | density of brine inclusions | Eqn. 6
$\rho_{f}$ | kg/m3 | final ambient fluid density | directly measured
$S_{f}$ | g/kg | final ambient fluid salinity | Eqn. 6; needs $\rho_{f}$
$S_{a}$ | g/kg | salinity of the ambient fluid | Eqn. 13; needs $S_{o}$, $S_{f}$, $h$, $h_{o}$, $h_{f}$
$T_{a}$ |  ∘C | temperature of the ambient fluid | directly measured and logged
$\rho_{a}$ | kg/m3 | density of the ambient fluid | Eqn. 6; needs $S_{a}$

The mass of the fluid expelled into the beaker, $m_{e}$, was logged in time
for the duration of each experiment. The salinity of the expelled fluid
$S_{e}$ was calculated from direct measurements of its density with a
precision Anton-Paar densimeter and the equation of state for NaCl and fresh
tapwater given by ?,

$\rho=998.43+0.69722S+2.5201\times 10^{-4}S^{2},$ (6)

which was intended for salinities up to $S=260$ g/kg and temperatures between
0∘C and -20∘C from the data of ?; note that for this temperature range, the
maximum error from neglecting the temperature dependence of the density is 2%
(?).

The volume of the tank $V_{t}$ is constant and determined from accurate
measurements of the tank geometry. The initial ambient fluid density
$\rho_{o}$ is directly measured with a precision Anton-Paar densimeter. The
initial ambient fluid salinity $S_{o}$ is calculated from its density with
Eqn. 6.

Following ?, we estimate the density of the solid ice $\rho_{ice}$ with the
equation of state given by ?,

$\rho_{ice}=916.8-0.1403T_{ice},$ (7)

where $T_{ice}$ is the solid ice temperature in ∘C. As we do not have direct
measurements of temperature within the solid ice, we approximate it as the
mid-point between the freezing plate temperature $T_{FP}$ and the ambient
fluid–ice interface temperature $T_{liq}$. The interface temperature $T_{liq}$
is assumed to be the liquidus temperature for water with the ambient salinity
$S$, which is given by ? as,

$T_{liq}(S)=-5.92\times 10^{-2}S-9.37\times 10^{-6}S^{2}-5.33\times
10^{-7}S^{3},$ (8)

for the range of $S$ up to 230 g/kg. Thus,

$T_{ice}(S)=\frac{T_{FP}+T_{liq}(S)}{2},$ (9)

which is used in Eqn. 7 to give us the density of solid ice $\rho_{ice}$.

The thickness of the mushy layer $h$ increases as the ice grows; this
thickness was estimated from the high-resolution images captured throughout
the experiments. The light intensity discontinuity across the mushy layer and
liquid ambient boundary was identified with image analysis software and used
to calculate the tank-average thickness of the mushy layer. This mushy layer
thickness was multiplied by the relevant geometry of the tank to obtain an
estimate of the mushy layer volume $V_{m}$. This simple approach assumes there
is no three-dimensional structure in the mushy layer growth, which in reality
is not the case. To account for three-dimensionality in the mushy layer, a
small geometrical volume correction is applied (described below). Note that
the initial mushy layer thickness $h_{o}$ is not necessarily zero; the first
photo of the experiment is typically 30 minutes after the freezing plate
temperature has been adjusted. The final mushy layer thickness $h_{f}$ is the
thickness of the mushy layer when it has reached an equilibrium, which was
defined as the time when the average mushy layer thickness has not increased
by more than 0.5% relative to the previous image.

The salinity of the brine inclusions $S_{b}$ is assumed to be the salinity
corresponding to the boundary between mushy layer and aqueous solution for
that particular ice temperature, which is given by the liquidus curve. That
is, the equation for the liquidus (Eqn. 8) can be inverted to give brine
salinity as a function of ice temperature (∘C),

$S_{b}=-17.6T_{ice}-0.389T_{ice}^{2}-3.62\times 10^{-3}T_{ice}^{3}.$ (10)

With the brine salinity we can then estimate the brine density with Eqn. 6.
Note that this equation for brine salinity as a function of ice interior
temperature is specific for aqueous solutions of NaCl and water, and differs
from that employed by sea ice models (?),

$S_{b}=\left({10^{-3}-\frac{0.054}{T_{ice}}}\right)^{-1},$ (11)

which is specific for brine salinity in sea ice formed from actual seawater.

The final salinity of the ambient fluid $S_{f}$ is calculated from the final
density of the ambient fluid $\rho_{f}$, which is measured directly with a
precision Anton-Paar densimeter, and Equation 6. The evolving ambient fluid
density $\rho_{a}$, however, is not measured during the experiment, so the
evolving salinity of the ambient fluid $S_{a}$ needs to be estimated using
measurements of the initial and final ambient fluid salinities, $S_{o}$ and
$S_{f}$, respectively. The salinity of the ambient fluid evolves due to the
growth of the mushy layer and the brine rejection from the mushy layer; these
processes are not constant in time, so approximating the ambient fluid
salinity with a linear evolution between the initial and final salinities is
not valid. Thus, we follow the approach of ?, and relate the evolution of the
ambient salinity to the evolution of the mushy layer thickness $h$ by way of a
scaling factor; that is,

$H_{c}=\frac{h(t)-h_{o}}{h_{f}-h_{o}},$ (12)

where $h(t)$ is the mushy layer thickness at time $t$, and $h_{o}$ and $h_{f}$
are the initial and final mushy layer thicknesses. The estimated salinity of
the ambient fluid is then given by,

$S_{a}(t)=H_{c}\left(S_{f}-S_{o}\right)+S_{o}.$ (13)

The temperature of the ambient fluid $T_{a}$ is directly measured and logged
throughout the experiment. The evolving density of the ambient fluid
$\rho_{a}$ is then able to be calculated from the ambient fluid temperature
and salinity with an appropriate equation of state, for which we use Equation
6.

In summary, the mass and salt conversation methods require the knowledge of
several salinities, temperatures, densities, masses, volumes, and thicknesses.
Some are directly measured ($V_{t}$, $m_{e}$, $\rho_{e}$, $\rho_{o}$,
$\rho_{f}$, $T_{a}$), some are estimated from high-resolution photos ($h$,
$h_{o}$, $h_{f}$, $V_{m}$), and some are approximated from established and/or
adapted relationships ($S_{e}$, $S_{o}$, $\rho_{ice}$, $T_{ice}$, $T_{liq}$,
$S_{b}$, $\rho_{b}$, $S_{f}$, $S_{a}$, $\rho_{a}$).

##  4 Results

Figure 3: Photos of the equilibrated mushy layer for the top, side, and bottom
freezing plate orientations (left to right) for experiments with initial
ambient salinities of $S_{o}=33$ g/kg and freezing plate temperatures of
$T_{FP}=-20^{\circ}$C.

Once the experiments began, the mushy layer was observed to form on the
freezing plate within the first 30 mins. The initial distribution of the mushy
layer was uniform across the freezing plate, indicative of a uniform boundary
condition imposed by the heat exchanger. In time, the edges of the mushy layer
within $\sim$5 mm of the tank sidewalls and ice face became rounded,
suggesting a small amount of unintentional heat entering the tank from the
laboratory; for the bottom oriented case, this rounded mushy layer edge was
more obvious (Fig. 3). Also, when estimating the mushy layer volume from the
high-resolution photos, a minor geometric correction is applied to account for
the mushy layer having rounded edges. As the mushy layer continues to grow,
the ice face tended to develop low-mode three-dimensional structures in the
plane parallel to the freezing plate (vertical and into the page in Fig. 2).
The exact nature of these structures appears to depend on minor differences in
initial conditions since repeating experiments didn’t necessarily reproduce
identical shapes of mushy layer face structures. Experiments with bottom
oriented freezing plates exhibited relatively less three-dimensionality to the
mushy layer face, suggesting the structure is perhaps related to the gravity
drainage process in the top and side orientations. Nevertheless, the amplitude
of these features were small relative to the mushy layer thickness; as such,
the mushy layer thickness and volume estimates were developed by taking a tank
average of all ice edges visible. The detection limit of this method is given
by the thickness per pixel, which is approximately 0.7 mm/pixel. Considering
the complications associated with the rounded edges and three-dimensionality
corrections, generous uncertainties (up to 5%, depending on the extent of the
rounded edges and three-dimensionality) are applied to the mushy layer volume
estimates.

Figure 4: The time evolution of bulk mushy layer thickness for all
experiments; the columns indicate the different freezing plate temperatures
increasing from left to right, and the rows indicate the different freezing
plate orientations. The line colours represent the different initial ambient
salinities (see legend; in g/kg). The vertical dashes indicate the times that
photos were taken, and their extent is indicative of the measurement
uncertainty. The circles represent the bulk mushy layer thickness when the
experiment has reached equilibrium.

Analysis of the high resolution photos allows us to quantify the bulk mushy
layer thickness and how it changes over the course of an experiment (Fig. 4).
For all experiments, the mushy layer thickness monotonically increases in
time, with the rate of thickness change tending to decrease in time. The
responses of the mushy layer thickness to the initial ambient salinity and
freezing plate temperatures are intuitive and consistent for the top and side
orientations; experiments with higher initial salinities and warmer freezing
plate temperatures exhibit reduced mushy layer thicknesses and growth rates.
The bottom orientation experiments exhibit less sensitivity to the initial
ambient salinity and freezing plate temperature, and faster growth rates.

Figure 5: The initial (hollow circles) and final (solid circles) ice bulk
salinity and interior temperature for the top, side, and bottom freezing plate
orientation (left to right). The horizontal lines through the solid circles
represent the range of final ice bulk salinities based on the two methods to
calculate the solid mass fraction. The yellow, white, and blue regions
represent the aqueous solution, mushy layer, and compositional mushy layer
regimes, respectively, as per Figure 1.

Figure 5 shows the evolution of the interior ice conditions. The ice interior
temperature characterises the bulk thermal conditions of the ice. The ice bulk
salinity reflects the ratio of the total mass of salt contained in brine
within the mushy layer (in grams) to the total mass of the mushy layer and
brine mixture (in kilograms); this is the salinity of the aqueous solution
that results from the melting of the mushy layer. Here, the initial ice bulk
salinity is assumed to be that of the initial ambient fluid salinity $S_{o}$,
and the initial ice interior temperature is that of the freezing plate
$T_{FP}$. The final ice interior temperature $T^{\star}$ is given by Eqn. 9,
where we use the freezing plate temperature $T_{FP}$, and the liquidus
temperature $T_{liq}$ based on the final ambient salinity $S_{f}$. The final
ice bulk salinity $S^{\star}$ is calculated using the solid mass fraction
$\phi$ and the brine salinity $S_{b}$, which is given by Eqn. 10 using the
final ice interior temperature $T^{\star}$; the relationship between the ice
bulk salinity, brine salinity, and solid mass fraction is,

$\phi=1-\frac{S^{\star}}{S_{b}}\quad\implies\quad
S^{\star}=\left(1-\phi\right)S_{b}.$ (14)

Here for $\phi$ we use the average value of $\phi_{M}$ and $\phi_{S}$, and
provide an indication of the differences of ice bulk salinities estimated by
the two different methods (horizontal lines through the final ice bulk
salinities).

The final ice interior temperature is, by definition, always warmer than the
initial ice interior temperature. The final ice bulk salinity is fresher that
the initial ice bulk salinity for the top and side orientations; for the
bottom orientation, however, the initial and final ice bulk salinities are
virtually the same. This freshening of the ice bulk salinity reflects the
gravity drainage of salty brine from the mushy layer in the top and side
orientations, which is a process that is not able to occur in the bottom
orientation experiments.

Figure 6: The time evolution of the solid mass fraction for all experiments;
the columns indicate the different freezing plate temperatures increasing from
left to right, and the rows indicate the different freezing plate
orientations. The line colours represent the different initial ambient
salinities (see legend; in g/kg). The upright triangles are the solid mass
fractions calculated with the salt conservation method $\phi_{S}$; the
inverted triangles are the solid mass fractions calculated with the mass
conservation method $\phi_{M}$; the lines follow the average value of the two.
The circles represent the solid mass fraction when the experiment has reached
equilibrium. The horizontal lines indicate the predicted solid mass fraction
$\phi$ based on the initial ice bulk salinity (approximated by $S_{o}$) and
interior ice temperature $T_{ice}$ (Eqns. 10 & 14).

The solid mass fractions calculated by the mass and salt conservation methods
over the course of each experiment are shown in Figure 6. In general, the
solid mass fractions decrease for increasing ambient salinities, and warming
freezing plate temperatures. The experiments with top and side orientations
exhibit solid mass fractions that increase in time, and the rate at which
these solid mass fractions increase tends to increase with ambient salinity.
The bottom orientation experiments tend to remain near their initial solid
mass fractions.

The initial ice bulk salinity (approximated by the initial ambient fluid
salinity $S_{o}$) and the ice interior temperature $T_{ice}$ can be used with
Equations 10 & 14 to calculate a predicted solid mass fraction, which are
shown in Figure 6 as horizontal lines. These predicted solid mass fractions
exhibit good agreement with the measured $\phi$ during early stages, and
throughout for the bottom orientation experiments. Interestingly, while the
top and bottom orientation experiments tend to have larger $\phi_{S}$ relative
to $\phi_{M}$, the side orientation experiments have larger $\phi_{M}$; the
reason for this is unknown.

Figure 7: Measured solid mass fractions at equilibrium plotted by their
respective freezing plate temperatures $T_{FP}$ and initial ambient salinities
$S_{o}$ (top row), and their respective final ice interior temperatures
$T^{\star}$ and bulk salinities $S^{\star}$ (bottom row). The background
colourmap indicates the solid mass fractions predicted by Equations 8, 9, 10 &
14 using a range ice bulk salinities and interior temperatures. The contours
are at $\phi=0.1$ intervals, and the magenta line represents the liquidus
curve.

The equilibrated solid mass fractions are sensitive to the salinity and
temperature conditions of the system. Figure 7 shows the measured solid mass
fractions of the experiments plotted by their respective initial ambient
salinities $S_{o}$ and freezing plate temperatures $T_{FP}$ (top row), and
their final ice bulk salinities $S^{\star}$ and interior temperatures
$T^{\star}$ (bottom row). These plots include the solid mass fractions
predicted by Equations 8, 9, 10 & 14 for a range of ambient salinities and
freezing plate temperatures (top row) and ice bulk salinities and /interior
temperatures (bottom row).

In general, there is good agreement between the measured and predicted solid
mass fractions; they exhibit consistent behaviours in terms of their relative
sensitivities to the salinity and temperature conditions. The agreement is
improved when using the final ice bulk salinities and interior temperatures
instead of the initial ambient salinities and freezing plate temperatures. The
agreement between the measured and predicted solid mass fractions is best for
experiments with smaller ice bulk salinities; this makes good physical sense
as the solid mass fraction converges to $\phi=1$ as system freshens. The
experiments with bottom freezing plate orientations exhibit the best agreement
for the initial conditions; this reflects the fact that brine drainage does
not occur in this configuration, such that the equilibrium ice conditions are
well represented by the known initial conditions. We hypothesise that the
leading sources of the minor disagreement between the measured and predicted
is uncertainty in regards to the final ice interior temperature (which we
assume to be the mid-point between the freezing plate and liquidus
temperatures; Eqn. 9), and variations throughout the mushy layer that are not
well represented by the bulk approach employed here.

The sensitivities of the measured solid mass fractions with respect to the
final ice bulk salinities and ice interior temperatures are shown in Figure 8,
plotted by the freezing plate temperature and initial ambient salinity,
respectively. These are calculated by taking the mean gradient of the measured
$\phi$ with respect to bulk ice salinity at the three different freezing plate
temperatures ($\partial\phi/\partial{S^{\star}}|_{T_{FP}}$; left panel), and
the mean gradient of the measured $\phi$ with respect to ice interior
temperature at the five different initial ambient salinities
($\partial\phi/\partial{T^{\star}}|_{S_{o}}$; right panel). Also included are
the predicted sensitivities derived from the solid mass fraction predictions
shown in Figure 7; for these we use values of
$\partial\phi/\partial{S^{\star}}|_{T_{FP}}$ at $S_{o}=33$ g/kg for the three
freezing plate temperatures, and the mean values of
$\partial\phi/\partial{T^{\star}}|_{S_{o}}$ at
$T_{FP}=\left(-20,-15,-10\right)^{\circ}$C for the five different initial
ambient salinities.

The solid mass fraction sensitivity to ice bulk salinity has a strong
dependence on the freezing plate temperature; for a given freezing plate
temperature, the rate of change of solid mass fraction with respect to ice
bulk salinity nearly doubles from
$\partial\phi/\partial{S^{\star}}|_{T_{FP}}\approx-0.005$ (g/kg)-1 at
$T_{FP}=-20^{\circ}$C to more than
$\partial\phi/\partial{S^{\star}}|_{T_{FP}}\approx-0.009$ (g/kg)-1 at
$T_{FP}=-10^{\circ}$C. This is reflected in the convergence of the predicted
solid mass fraction contours for warming freezing plate temperatures in Figure
7; this convergence occurs because $\phi$ goes from $\phi=1$ at zero
salinities to $\phi=0$ at the liquidus salinity, which decreases for
increasing temperatures. The sensitivity of the solid mass fraction to the ice
interior temperature depends on the initial ambient salinity; the rate of
change of solid mass fraction with ice interior temperature increases by a
factor of 5 between initial ambient salinities of $S_{o}\approx 18$ g/kg to
$S_{o}\approx 102$ g/kg. There is good agreement between the measured
sensitivities and those predicted by the relationship in Equations 8, 9, 10 &
14.

The sensitivity analysis allows us to determine whether the solid mass
fraction is more responsive to the temperature or salinity properties of the
system. For the oceanographically realistic experiments with $S_{o}\approx 33$
g/kg, the sensitivity of solid mass fraction to ice interior temperature is
approximately
$\partial\phi/\partial{T^{\star}}|_{S_{o}}\approx-0.02^{\circ}$C-1, or -2.5%
per 1∘C increase. To achieve an equivalent magnitude change in the solid mass
fraction from a change in the ice bulk salinity requires a salinity increase
of approximately 4 g/kg at $T_{FP}=-20^{\circ}$C, or approximately 2 g/kg at
$T_{FP}=-10^{\circ}$C. When considering the relative ranges of temperature and
salinity variabilities at high latitudes, this sensitivity analysis suggests
that the solid mass fraction is temperature dominated. It follows that
increases in ice interior temperature will lead to a reduction in the ice
solid mass fraction.

Figure 8: The sensitivities of measured solid mass fraction to the ice bulk
salinity $\partial\phi/\partial{S^{\star}}|_{T_{FP}}$ (left) and the ice
interior temperature $\partial\phi/\partial{T^{\star}}|_{S_{o}}$ (right) for
given freezing plate temperatures and initial ambient salinities,
respectively. The predicted sensitivities given by the relationships in
Equations 8, 9, 10 & 14 are also included (blue stars).

##  5 Discussion & Conclusions

The experiments demonstrate how the mushy layer growth rate and solid mass
fraction are affected by the freezing plate temperature, ambient salinity, and
freezing orientation. Our approach has allowed for a thorough exploration of a
wide range of thermodynamic conditions and freezing orientations with a single
apparatus and common methodology. The results exhibit good agreement with
existing theories for how the solid mass fraction depends on ice bulk salinity
and interior temperature (e.g., ???). Mushy layer growth rate is largest for
cooler freezing plate temperatures, fresher ambient salinities, and bottom
oriented freezing configurations. The equilibrated mushy layer solid mass
fractions are largest for cooler freezing plate temperatures and fresher
ambient salinities, and smallest for the bottom oriented freezing
configuration where brine drainage is unable to occur. The sensitivity of the
solid mass fraction to the ice temperature and salinity is explored, and it is
found that for oceanographically-realistic conditions the solid mass fraction
is most sensitive to changes in temperature; every 1∘C increase in freezing
temperatures will result in a 2.5% decrease of solid mass fraction.

The bottom oriented freezing experiments exhibit substantially different mushy
layer behaviour and solid mass fraction evolution, which results from the
absence of brine drainage. The top and side oriented freezing experiments
exhibit very little difference; they have similar solid mass fraction
evolutions, final values, and sensitivities for the oceanographically relevant
range of salinities (i.e., $S_{o}=\left(18,33\right)$ g/kg). This finding
suggests that the introduction of floe size distribution capabilities in
modern numerical sea ice models (e.g., CICE6 v6.1.2), and the associated
increase in lateral sea ice dynamics, should not require distinct
parameterisations for solid mass fractions arising from vertical versus
lateral sea ice growth.

The implications of the relatively high thermal sensitivity of solid mass
fraction are worthy of closer attention. Increases in the freezing
temperatures of sea ice, which are observed and predicted to continue under
realistic climate forcing scenarios (e.g., ?), will reduce the solid mass
fraction of sea ice. This solid mass fraction reduction will decrease the sea
ice brine rejection and associated salt flux into the ocean, and reduce the
capacity of the sea ice field to freshen the high latitude oceans during melt.
Decreasing the brine rejection salt flux will ultimately lead to a freshening
of the high-latitude waters that feed the deep, bottom and abyssal waters of
the global oceans, thereby having the potential to reduce the strength of the
global thermohaline overturning circulation. A reduced solid mass fraction
equates to an increased sea ice bulk salinity; when this saltier sea ice
melts, the surrounding surface waters are freshened less than they would
otherwise be with higher solid mass fraction sea ice. This increase in the
high-latitude surface salinity would tend to weaken the surface stratification
of these salt stratified waters (i.e., beta-oceans; ?), and lead to a reduced
sea ice growth rate. All of these implications share a common consequence of
tending to reduce the role of sea ice in the climate system. Investigating the
potential feedbacks caused by solid mass fraction variations in coupled
ocean–sea ice and climate models is a primary topic for future research.

##  6 Acknowledgments

We wish to thank Angus Rummery and Tony Beasley for the construction of the
apparatus and laboratory assistance.

## Bibliography

  * Anderson D and Guba P (2020) Convective phenomena in mushy layers. _Annual Review of Fluid Mechanics_ , 52, 93–119 (10.1146/annurev-fluid-010719-060332)
  * Budyko M (1969) The effect of solar radiation variations on the climate of the earth. _Tellus_ , 21, 611–619
  * Eayrs C, Holland D, Francis D, Wagner T, Kumar R and Li X (2019) Understanding the seasonal cycle of antarctic sea ice extent in the context of longer-term variability. _Reviews of Geophysics_ , 57, 1037–1064 (10.1029/2018RG000631)
  * Feltham D, Untersteiner N, Wettlaufer J and Worster M (2006) Sea ice is a mushy layer. _Geophysical Research Letters_ , 33, L14501 (10.1029/2006GL026290)
  * Hunke E, Lipscomb W, Turner A, Jeffery N and Elliot S (2015) Cice: the los alamos sea ice model documentation and software user’s manual version 5.1. _Los Alamos National Laboratory_
  * Huppert H (1990) The fluid mechanics of solidification. _Journal of Fluid Mechanics_ , 212, 2090240 (10.1017/S0022112090001938)
  * Kraitzman N, Promislow K and Wetton B (2022) Slow migration of brine inclusions in first-year sea ice. _Society for Industrial and Applied Mathematics: Journal on Applied Mathematics_ , 82, 1470–1494 (10.1137/21M1440244)
  * Manabe S and Stouffer R (1979) A co2-climate sensitivity study with a mathematical model of the global climate. _Nature_ , 282, 491–493
  * Notz D (2005) _Thermodynamic and fluid-dynamical processes in sea ice (Ph.D. Thesis)_. University of Cambridge, Cambridge, United Kingdom
  * Notz D and Worster M (2009) Desalination process of sea ice revisited. _Journal of Geophysical Research: Oceans_ , 114, C05006 (10.1029/2008JC004885)
  * Peppin S, Aussillous P, Huppert H and Worster M (2007) Steady-state mushy layers: Experiments and theory. _Journal of Fluid Mechanics_ , 570, 69–77 (10.1017/S0022112006003028)
  * Post E, Alley R, Christensen T, Macias-Fauria M, Forbes B, Gooseff M, Iler A, Kerby J, Laidre K, Mann M, Olofsson J, Stroeve J, Ulmer F, Virginia R and Wang M (2019) The polar regions in a 2∘c warmer world. _Science Advances_ , 5, eaaw9883 (10.1126/sciadv.aaw9883)
  * Pounder E (1965) _The Physics of Ice_. Pergamon Press, Oxford, United Kingdom
  * Purich A and Doddridge E (2023) Record low antarctic sea ice coverage indicates a new sea ice state. _Nature: Communications Earth & Environment_, 4, 314 (10.1038/s43247-023-00961-9)
  * Scotti K, Kearney L, Burns J, Ocana M, Duros L, Shelhamer A and Dunard D (2019) The effect of solidification direction with respect to gravity on ice-templated tio2 microstructures. _Journal of the European Ceramic Society_ , 39, 3180–3193 (10.1016/j.jeurceramsoc.2019.04.007)
  * Stewart K and Haine T (2016) Thermobaricity in the transition zones between alpha and beta oceans. _Journal of Physical Oceanography_ , 46, 1805–1821 (10.1175/JPO-D-16-0017.1)
  * Vancoppenolle M, Madec G, Thomas M and McDougall T (2019) Thermodynamics of sea ice phase composition revisited. _Journal of Geophysical Research: Oceans_ , 124, 615–634 (10.1029/2018JC014611)
  * Weast R (1971) _Handbook of Chemistry and Physics_. Chemical Rubber Company, Cleveland, Ohio, USA
  * Weeks W (2010) _On Sea Ice_. University of Alaska Press, Fairbanks, Alaska
  * Weeks W and Ackley S (1986) _The Growth, Structure, and Properties of Sea Ice. In: The Geophysics of Sea Ice, N. Untersteiner (ed.)_ , 9–164. Plenum Press, New York, New York
  * Wells A, Hitchen J and Parkinson J (2019) Mushy-layer growth and convection, with application to sea ice. _Philosophical Transactions of the Royal Society A_ , 377, 20180165 (10.1098/rsta.2018.0165)
  * Wettlaufer J, Worster M and Huppert H (1997) Natural convection during solidification of an alloy from above with application to the evolution of sea ice. _Journal of Fluid Mechanics_ , 237, 649–669 (10.1017/S0022112092003562)
  * Worster M (1997) Convection in mushy layers. _Annual Reviews of Fluid Mechanics_ , 29, 910122 (10.1146/annurev.fluid.29.1.91)
  * Worster M and Rees Jones D (2015) Sea ice thermodynamics and brine drainage. _Philosophical Transactions of the Royal Society A_ , 373, 20140166 (10.1098/rsta.2014.0166)