-
Notifications
You must be signed in to change notification settings - Fork 9
/
eqice.tex
825 lines (784 loc) · 34.3 KB
/
eqice.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
\section{Ice Model Formulation}
\label{Iphys}
The sea-ice component of ROMS is a combination of the
elastic-viscous-plastic (EVP) rheology
\citep{Hunke97,Hunke_2001} and simple one-layer
ice and snow thermodynamics with a molecular sublayer under the ice
\citep{Mellor89}. It is tightly coupled, having
the same grid (Arakawa-C) and timestep as the ocean and sharing the same
parallel coding structure for use with MPI or OpenMP \citep{Budgell05}.
\subsection{Dynamics}
The momentum equations describe the change in ice/snow velocity due
to the combined effects of the Coriolis force, surface ocean tilt,
air and water stress, and internal ice stress:
(\ref{ist1}) and (\ref{ist2}):
\begin{align}
\frac{d}{ dt}(A \rho_i h_i u)
& = A \rho_i h_i fv - A \rho_i h_i g \frac{\partial \zeta_w }{ \partial x} +
A (\tau_a^x + \tau_w^x) + {\cal F}_x
\label{ist1} \\
\frac{d}{ d t}(A \rho_i h_i v)
& = - A \rho_i h_i fu - A \rho_i h_i g \frac{\partial \zeta_w }{ \partial y} +
A (\tau_a^y + \tau_w^y) + {\cal F}_y.
\label{ist2}
\end{align}
In this model, we neglect the nonlinear advection terms as well as
the curvilinear terms in the internal ice stress.
Nonlinear formulae are used for both the ocean-ice and air-ice surface
stress:
\begin{align}
\vec{\tau}_a & = \rho_a C_a | \vec{V}_{10} | \vec{V}_{10} \\
C_a & = \frac{1 }{ 2} C_d \left[ 1 - \cos( 2 \pi \min(h_i+.1, .5)
\right] \\
\vec{\tau}_w & = \rho_w C_w | \vec{v}_w - \vec{v} |
( \vec{v}_w - \vec{v}) .
\end{align}
The force due to the internal ice stress is given by the divergence of
the stress tensor $\sigma$. The rheology is given by the stress-strain
relation of the medium. We would like to emulate the viscous-plastic
rheology of \citet{Hibler79}:
\begin{equation}
\sigma_{ij} = 2 \eta \dot \epsilon_{ij} + (\zeta - \eta) \dot
\epsilon_{kk} \delta_{ij} - \frac{P }{ 2} \delta_{ij}
\end{equation}
\begin{equation}
\dot \epsilon_{ij} \equiv \frac{1 }{ 2} \left( \frac{\partial u_i
}{
\partial x_j} + \frac{\partial u_j }{ \partial x_i} \right)
\end{equation}
where the nonlinear viscosities are given by:
\begin{equation}
\zeta = \frac{ P }{ 2 \left[ (\epsilon^2_{11} +
\epsilon^2_{22} ) ( 1 + 1/e^2 ) + 4 e^{-2} \epsilon^2_{12}
+ 2 \epsilon_{11} \epsilon_{22} ( 1 - 1/e^2 ) \right] ^{1/2} }
\end{equation}
\begin{equation}
\eta = \frac{ \zeta }{ e^2 }.
\end{equation}
Hibler's ice strength is given by:
\begin{equation}
P = P^* A h_i e^{-C(1-A)}
\end{equation}
while \citet{Overland_1988} advocate a nonlinear strength:
\begin{equation}
P = P^* A h_i^2 e^{-C(1-A)}
\end{equation}
in which $P^*$ now depends on the grid spacing. Both options are available
in ROMS.
We would also like to have an explicit model that can be solved
efficiently on parallel computers. The EVP rheology has a tunable
coefficient $E$ (the Young's modulus) which can be chosen to make
the elastic term small compared to the other terms. We rearrange the VP
rheology:
\begin{equation}
\frac{1 }{ 2 \eta} \sigma_{ij} + \frac{\eta - \zeta }{ 4 \eta \zeta}
\sigma_{kk} \delta_{ij} + \frac{P }{ 4 \zeta} \delta_{ij} = \dot
\epsilon_{ij}
\end{equation}
then add the elastic term:
\begin{equation}
\frac{1 }{ E} \frac{\partial \sigma_{ij} }{ \partial t} +
\frac{1 }{ 2
\eta} \sigma_{ij} + \frac{\eta - \zeta }{ 4 \eta \zeta} \sigma_{kk}
\delta_{ij} + \frac{P }{ 4 \zeta} \delta_{ij} = \dot \epsilon_{ij}
\end{equation}
Much like the ocean model, the ice model has a split timestep. The
internal ice stress term is updated on a shorter timestep so as to
allow the elastic wave velocity to be resolved.
Once the new ice velocities are computed, the ice tracers can be
advected using the MPDATA scheme \citep{Smolark90}. The tracers in
this case are the ice thickness, ice concentration, snow thickness,
internal ice temperature, and surface melt ponds. The continuity
equations describing the evolution of these parameters (equations
(\ref{ist3a})--(\ref{ist3b})) also include thermodynamic terms ($S_h$,
$S_s$ and $S_A$), which will be described in \S\ref{Growth}:
\begin{align}
\frac{\partial A h_i }{ \partial t} & =
- \frac{\partial (u A h_i) }{ \partial x} -
\frac{\partial (v A h_i) }{ \partial y}
+ S_h + {\cal D}_h
\label{ist3a} \\
\frac{\partial A h_s }{ \partial t} & =
- \frac{\partial (u A h_s) }{ \partial x} -
\frac{\partial (v A h_s) }{ \partial y}
+ S_s + {\cal D}_s
\label{ist3c} \\
\frac{\partial A }{ \partial t} & =
- \frac{\partial (uA) }{ \partial x} - \frac{\partial (vA) }{ \partial y}
+ S_A + {\cal D}_A \qquad \qquad 0 \leq A \leq 1 .
\label{ist3b}
\end{align}
The first two equations represent the conservation of ice and snow.
Equation (\ref{ist3b}) is discussed in some detail in \citet{Mellor89}, but
represents the advection of ice blocks in which no ridging occurs as
long as there is any open water.
%An optional ridging term can be added
%\citep{Gray96}:
%\begin{equation}
% \frac{\partial A }{ \partial t} =
% - \frac{\partial (uA) }{ \partial x} - \frac{\partial (vA) }{ \partial y}
% - A \alpha(A) \, \nabla \cdot \vec{v} \, H(-\nabla \cdot \vec{v})
% + S_A + {\cal D}_A \qquad \qquad 0 \leq A \leq 1 .
%\end{equation}
%where $\alpha(A)$ is an arbitrary function such that $\alpha(0) = 0$,
%$\alpha(1) = 1$, and $0 \leq \alpha(A) \leq 1$. The ridging term leads
%to an increase in $h_i$ under convergent flow as would be produced by
%ridging. The function $\alpha(A)$ should be chosen so that it is near
%zero until the ice concentration is large enough that ridging is
%expected to occur, then should increase smoothly to one.
%
The symbols used in these equations along with the values for the
constants are listed in Table \ref{icemomvars}.
\begin{table}[pt]
\hspace{9.5 mm}
\vbox{
\begin{tabular}{|c|c|l|} \hline
Variable & Value & Description \\ \hline
$A(x,y,t)$ && ice concentration \\
% $\alpha(A)$ && ridging function \\
$C_a$ && nonlinear air drag coefficient \\
$C_d$ & $2.2 \times 10^{-3}$ & air drag coefficient \\
$C_w$ & $10 \times 10^{-3}$ & water drag coefficient \\
(${\cal D}_h, {\cal D}_s, {\cal D}_A$) && diffusion terms \\
$\delta_{ij}$ && Kronecker delta function \\
$E$ && Young's modulus \\
$e$ & 2 & eccentricity of the elliptical yield curve \\
$\epsilon_{ij}(x,y,t)$ && strain rate tensor \\
$\eta(x,y,t)$ && nonlinear shear viscosity \\
$f(x,y)$ && Coriolis parameter \\
(${\cal F}_x, {\cal F}_y$) && internal ice stress \\
$g$ & $9.8\, m\, s^{-2}$ & acceleration of gravity \\
$H$ && Heaviside function \\
$h_i(x,y,t)$ && ice thickness of ice-covered fraction \\
$h_o$ & 1 m & ice cutoff thickness \\
$h_s(x,y,t)$ && snow thickness on ice-covered fraction \\
$M(x,y,t)$ && ice mass (density times thickness) \\
$P(x,y,t)$ && ice pressure or strength \\
($P^*, C$) & ($2.75 \times 10^4, 20$) & ice strength parameters \\
($S_h, S_s, S_A$) && thermodynamic terms \\
$\sigma_{ij}(x,y,t)$ && stress tensor \\
$\vec{\tau}_a$ && air stress \\
$\vec{\tau}_w$ && water stress \\
($u,v$) && the ($x,y$) components of ice velocity $\vec{v}$ \\
($\vec{V}_{10}, \vec{v}_w$) &&
10 meter air and surface water velocities \\
($\rho_a, \rho_w$) & ($1.3\, kg\, m^{-3}, 1025\, kg\, m^{-3}$) & air
and water densities \\
$\zeta(x,y,t)$ && nonlinear bulk viscosity \\
$\zeta_w(x,y,t)$ && height of the ocean surface \\
\hline
\end{tabular}
}
\caption{Variables used in the ice momentum equations.}
\label{icemomvars}
\end{table}
Note that Hibler's $h_I$ variable is equivalent to our $A h_i$
combination - his $h_I$ is the average thickness over the whole
gridbox while our $h_i$ is the average thickness over the ice-covered
fraction of the gridbox.
\subsubsection{Landfast ice}
\label{Landfast_ice}
The Arctic ocean has many shallow shelves, on which landfast ice
can form every winter. The model as described thus far has no way to
reproduce the landfast ice, but \citet{Lemieux_2015} describes a
landfast ice parameterization. Basal stress terms $\tau_b$ are added to
equations (\ref{ist1}) and (\ref{ist2}) representing a bottom drag
applied to the deepest ice keels. Included is a parameterization of the
ice thickness distribution so that even a single ice-category model can
use the landfast ice option. Figure \ref{f_landfast} shows a cartoon of
an ice floe with a keel.
\begin{figure}[ht]
\begin{center}
\definecolor{cc87137}{RGB}{200,113,55}
\definecolor{c2ad4ff}{RGB}{42,212,255}
\definecolor{c55ddff}{RGB}{85,221,255}
\definecolor{c0072ff}{RGB}{0,114,255}
\begin{tikzpicture}[y=0.80pt, x=0.80pt, yscale=-1.000000, xscale=1.000000, inner sep=0pt, outer sep=0pt]
\begin{scope}[shift={(0,-485.43307)}]
\path[fill=cc87137,rounded corners=0.0000cm] (149.3704,752.5492) rectangle
(449.3704,772.5492);
\path[draw=c2ad4ff,fill=c55ddff,line join=miter,line cap=butt,even odd rule,line
width=0.800pt] (150.0000,602.3622) -- (150.0000,622.3622) --
(250.0000,622.3622) -- (260.0000,741.9196) -- (270.0000,622.3622) --
(360.0000,622.3622) -- (360.0000,602.3622) -- (270.0000,602.3622) --
(260.0000,592.3622) -- (250.0000,602.3622) -- cycle;
\path[draw=c0072ff,line join=miter,line cap=butt,miter limit=4.00,even odd
rule,line width=1.600pt] (360.0000,606.7695) .. controls (365.0000,606.7695)
and (374.1050,600.8756) .. (377.6873,601.0614) .. controls (381.2696,601.2471)
and (379.9588,603.7348) .. (395.6076,607.9057) .. controls (406.2484,608.4328)
and (407.0792,603.2019) .. (417.7482,602.1621) .. controls (430.6436,599.4896)
and (427.8328,609.1786) .. (441.1189,609.7014) .. controls (452.2796,609.3679)
and (447.9592,604.7616) .. (458.6417,600.6002) .. controls (461.0764,598.8109)
and (473.4825,605.7417) .. (472.3709,604.8806);
\path[draw=c0072ff,line join=miter,line cap=butt,miter limit=4.00,even odd
rule,line width=1.600pt] (150.0000,604.2510) .. controls (148.3100,605.5942)
and (142.4698,606.3718) .. (140.0415,606.6003) .. controls (125.5603,609.5491)
and (124.8491,595.8552) .. (116.5882,605.9987) .. controls (111.6403,609.6048)
and (104.8953,602.7983) .. (99.3704,605.4365);
\end{scope}
\path[draw=black,dash pattern=on 2.40pt off 0.80pt,line join=miter,line
cap=butt,miter limit=4.00,even odd rule,line width=0.800pt] (150.0000,35.4486)
-- (150.0000,95.5276);
\path[draw=black,dash pattern=on 2.40pt off 0.80pt,line join=miter,line
cap=butt,miter limit=4.00,even odd rule,line width=0.800pt] (450.0000,36.7078)
-- (450.0000,98.0460);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(280.8509,47.3375) -- (150.6296,47.3375);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(280.2213,76.7816) -- (150.6296,76.7816);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(296.0000,77.4112) -- (360.5953,77.4112);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(308.0000,47.3375) -- (450.0000,47.3375);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(280.0000,186.9291) -- (280.0000,136.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(280.0000,206.9291) -- (280.0000,256.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(210.0000,86.9291) -- (210.0000,116.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(210.0000,166.9291) -- (210.0000,136.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(250.0000,86.9291) -- (250.0000,106.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(250.0000,146.9291) -- (250.0000,116.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(360.0000,86.9291) -- (361.2592,120.0034);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(360.0000,166.9291) -- (360.0000,136.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(149.6707,180.0034) -- (149.6707,120.5593);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(149.7818,197.4850) -- (150.4114,265.7437);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(210.0000,86.9291) -- (210.0000,116.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(210.0000,166.9291) -- (210.0000,136.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(250.0000,86.9291) -- (250.0000,106.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(250.0000,146.9291) -- (250.0000,116.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(360.0000,86.9291) -- (361.2592,120.0034);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(360.0000,166.9291) -- (360.0000,136.9291);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(219.2597,121.9665) -- (249.2597,121.9665);
\path[->,>=stealth,draw=black,line join=miter,line cap=butt,even odd rule,line width=0.800pt]
(302.1470,121.8554) -- (272.1470,121.8554);
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(269.3704,202.5219) node[above right] (text7731) {$h_{rb}$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(362.4447,135.0403) node[above right] (text7735) {$h_d$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(283.7776,52.5219) node[above right] (text7739) {$\Delta x$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(284.4073,80.6330) node[above right] (text7743) {$l_i$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(143.1481,194.4107) node[above right] (text7747) {$h_w$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(193.1481,108.8180) node[above right] (text7751) {$h_l$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(231.8888,113.1515) node[above right] (text7755) {$h_{rt}$};
\path[fill=black,line join=miter,line cap=butt,line width=0.800pt]
(254.3626,127.4611) node[above right] (text4378) {$\beta l_i$};
\end{tikzpicture}
\caption{Cartoon of an ice floe with a deep keel.}
\label{f_landfast}
\end{center}
\end{figure}
\begin{table}[pt]
\begin{center}
%\hspace{9.5 mm}
\vbox{
\begin{tabular}{|c|c|l|} \hline
Variable & Value & Description \\ \hline
$h_w$ && water depth \\
$h_{rt}$ && ridge height above level ice \\
$h_{rb}$ && draft of keels below level ice bottom \\
$h_d$ && draft of ice below sea surface \\
$l_i$ && length of ice floe within grid cell \\
$\Delta x$ && $x$-dimension of grid cell \\
$\beta$ & O(0.01) & ridged fraction of floe \\
$k_1$ & 8.0 & tunable parameter \\
$k_2$ & 15.0 & tunable parameter \\
$u_0$ & $5.0\times 10^{-5}$ & small velocity \\
$C_b$ & 20 & ice strength parameter \\
\hline
\end{tabular}
}
\end{center}
\caption{Variables used in the landfast ice parameterization.}
\label{fasticevars}
\end{table}
Introducing the tunable parameters $k_1$ and $k_2$, plus the
expression $h_c = A h_w/k_1$, the $u$-component of the basal
stress due to the dragging of ice keels is computed as:
\begin{equation}
\tau_{bu} =
\begin{cases}
0& \text{if $h_w \geq h_c$},\\
k_2 \left( \frac{u}{|\vec{v}| + u_0} \right)
(h_w-h_c)\exp^{-C_b(1-A)}& \text{if $h_w < h_c$}.
\end{cases}
\end{equation}
where $h_w$ is the water depth.
The expression for the $v$-component is similar.
\subsection{Thermodynamics}
\label{Growth}
The thermodynamics is based on calculating how much ice grows and
melts on each of the surface, bottom, and sides of the ice floes,
as well as frazil ice formation \citep{Mellor89}.
Once the ice tracers are advected, the ice concentration and
thickness are timestepped according to the terms on the right.
Equations (\ref{ist3a}) and (\ref{ist3b}) become:
\begin{align}
\frac{D A h_i }{D t}
& = \frac{\rho_o }{ \rho_i} \left[ A (W_{io} - W_{ai}) + (1-A) W_{ao}
+ W_{fr} \right]
\label{ist4a} \\
\frac{D A }{ D t}
& = \frac{\rho_o}{ \rho_i h_i} \left[ \Phi (1-A) W_{ao} + (1-A)
W_{fr} \right] \qquad \qquad 0 \leq A \leq 1 .
\label{ist4b}
\end{align}
The term $Ah_i$ is the ``effective thickness'', a measure of the ice
volume. Its evolution equation is simply quantifying the change in
the amount of ice. The ice concentration equation is more interesting in
that it provides the partitioning between ice melt/growth on the sides
vs. on the top and bottom. The parameter $\Phi$ controls this and has
differing values for ice melt and retreat. In principle, most of the ice
growth is assumed to happen at the base of the ice while rather more of
the melt happens on the sides of the ice due to warming of the water in
the leads.
The heat fluxes through the ice are based on a simple one layer
\citet{Semtner76a} type model with snow on top. The temperature is
assumed to be linear within the snow and within the ice. The ice contains
brine pockets for a total ice salinity of 3.2 or the surface salinity,
whichever is less. The surface ocean temperature and salinity is half a
$dz$ below the surface. The water right below the surface is assumed to
be at the freezing temperature; a logarithmic boundary layer is computed
having the temperature and salinity matched at freezing.
Here, the $W$ variables are the freeze or
melt rates as shown in Fig.\ \ref{fm+k} and Table \ref{thermvar}. The
frazil ice growth $W_{fr}$ will be discussed further in
\S\ref{frazil}---note that it contributes to changes in $A$ as well as
to changes in $h_i$. The other term that contributes to $A$ is
$W_{ao}$. This term includes a factor $\Phi$ which Mellor and Kantha
set to different values depending on whether ice is melting or
freezing:
\begin{align}
\Phi & = 4.0 \qquad\qquad W_{ao} \geq 0 \\
\Phi & = 0.5 \qquad\qquad W_{ao} < 0
\end{align}
\begin{figure}[ht]
\setlength{\unitlength}{0.00083300in}%
%
%\begin{picture}(6527,2507)(270,-5042)
\begin{picture}(6591,3076)(835,-5072)
%\thicklines
\put(1501,-2911){\line( 1, 0){3975}}
\put(1501,-4111){\line( 1, 0){3975}}
\put(1501,-2911){\line( 0,-1){1200}}
\put(6976,-4111){\line( 1, 0){1050}}
\put(6976,-2941){\line( 1, 0){1050}}
\put(6976,-2711){\line( 0,-1){1400}}
\put(6976,-2711){\line( 1, 0){1050}}
\put(5476,-2911){\line( 0,-1){1200}}
\put(2401,-4811){\vector( 0, 1){700}}
\put(3376,-4861){\vector( 0, 1){750}}
\put(6226,-3586){\vector( 0, 1){675}}
\put(4426,-4036){\vector( 0,-1){525}}
\put(2596,-3376){\vector( 1, 1){450}}
\put(7526,-2211){\vector( 0,-1){500}}
\put(7490,-3383){\vector( 1,2){200}}
\put(6121,-3786){\makebox(0,0)[lb]{$W_{ao}$}}
\put(4336,-4756){\makebox(0,0)[lb]{$W_{ro}$}}
\put(2416,-3586){\makebox(0,0)[lb]{$W_{ai}$}}
\put(3241,-5056){\makebox(0,0)[lb]{$W_{io}$}}
\put(2296,-5056){\makebox(0,0)[lb]{$W_{fr}$}}
\put(7400,-2180){\makebox(0,0)[lb]{$W_{s}$}}
\put(7276,-3611){\makebox(0,0)[lb]{$W_{sm}$}}
\end{picture}
%\begin{picture}(0,0)(5346,0)
\begin{picture}(0,0)(5961,-30)
\includegraphics{pics/therm_mk}%
\end{picture}%
\caption{Diagram of the different locations where ice melting and
freezing can occur.}
\label{fm+k}
\end{figure}
Similar to Eq. (\ref{ist4a}) is the snow equation:
\begin{equation}
{D A h_s \over D t} = \left[ A (W_{s} - W_{sm}) \right]
\end{equation}
where $W_s$ and $W_{sm}$ are the snowfall and snow melt rates,
respectively, in units of equivalent water. Currently all snow and
ice melt ends up in $W_{ro}$, though melt ponds could be added with
some effort.
\begin{table}
\hspace{9 mm}
\vbox{
\begin{tabular}{|c|c|l|} \hline
Variable & Value & Description \\ \hline
$\alpha_w$ & 0.10 & shortwave albedo of water \\
$\alpha_i$ & 0.60, 0.65 & shortwave albedo of wet, dry ice \\
$\alpha_s$ & 0.72, 0.85 & shortwave albedo of wet snow \\
$C_k$ && snow correction factor \\
$C_{pi}$ & 2093 J kg$^{-1}$ K$^{-1}$ & specific heat of ice \\
$C_{po}$ & 3990 J kg$^{-1}$ K$^{-1}$ & specific heat of water \\
$\epsilon_w$ & 0.97 & longwave emissivity of water \\
$\epsilon_i$ & 0.97 & longwave emissivity of ice \\
$\epsilon_s$ & 0.99 & longwave emissivity of snow \\
$E(T,r)$ && enthalpy of the ice/brine system \\
$F_T\!\uparrow$ && heat flux from the ocean into the ice \\
$h^\star$ && test for thick snow pushing ice surface under water \\
$H\!\downarrow$ && sensible heat \\
$I_o$ & part of $(1-\alpha_i)SW\!\!\downarrow$ & shortwave radiation entering ice \\
$i_w$ && fraction of the solar heating transmitted \\
&& through a lead into the water below \\
$k_i$ & 2.04 W m$^{-1}$ K${^-1}$ & thermal conductivity of ice \\
$k_s$ & 0.31 W m$^{-1}$ K${^-1}$ & thermal conductivity of snow \\
$L_i$ & 302 MJ m$^{-3}$ & latent heat of fusion of ice \\
$L_s$ & 110 MJ m$^{-3}$ & latent heat of fusion of snow \\
$LE\!\downarrow$ && latent heat \\
$LW\!\!\downarrow$ && incoming longwave radiation \\
$m$ & $-0.054^\circ$C & coefficient in linear $T_f(S) = mS$ equation \\
$\Phi$ && contribution to $A$ equation from freezing water \\
$Q_{ai}$ && heat flux out of the snow/ice surface \\
$Q_{ao}$ && heat flux out of the ocean surface \\
$Q_{i2}$ && heat flux up out of the ice \\
$Q_{io}$ && heat flux up into the ice \\
$Q_{s}$ && heat flux up through the snow \\
$r$ & $0 \le r \le 0.2 $ & brine fraction in ice \\
$\rho_i$ & 910 $m^3/kg$ & density of ice \\
$S_i$ & 3.2 & salinity of the ice \\
$SW\!\!\downarrow$ && incoming shortwave radiation \\
$\sigma$ & $5.67 \times 10^{-8}$ W m$^{-2}$ K$^{-4}$ &
Stefan-Boltzmann constant \\
$T_0$ && temperature of the bottom of the ice \\
$T_1$ && temperature of the interior of the ice \\
$T_2$ && temperature at the upper surface of the ice \\
$T_3$ && temperature at the upper surface of the snow \\
$T_f$ && freezing temperature \\
$T_{{\rm melt}\_i}$ & $mS_i$ & melting temperature of ice \\
$T_{{\rm melt}\_s}$ & 0$^\circ$ C & melting temperature of snow \\
$W_{ai}$ && melt rate on the upper ice/snow surface \\
$W_{ao}$ && freeze rate at the air/water interface \\
$W_{fr}$ && rate of frazil ice growth \\
$W_{io}$ && freeze rate at the ice/water interface \\
$W_{ro}$ && rate of run-off of surface melt water \\
$W_{s}$ && snowfall rate \\
$W_{sm}$ && snow melt rate \\
\hline
\end{tabular}
}
\caption{Variables used in the ice thermodynamics}
\label{thermvar}
\end{table}
\begin{figure}
\centerline{
\begin{picture}(0,0)%
\includegraphics{pics/therm_mk2}%
\end{picture}%
\setlength{\unitlength}{3947sp}%
%
\begin{picture}(6591,2608)(1468,-5057)
\put(6376,-3511){\makebox(0,0)[lb]{{$F_T$}}}
\put(3376,-4561){\makebox(0,0)[lb]{{$F_T$}}}
\put(3376,-4036){\makebox(0,0)[lb]{{$Q_{io}$}}}
\put(3376,-3436){\makebox(0,0)[lb]{{$Q_{i2}$}}}
\put(6376,-3061){\makebox(0,0)[lb]{{$Q_{ao}$}}}
\put(3376,-3136){\makebox(0,0)[lb]{{$Q_s$}}}
\put(3376,-2761){\makebox(0,0)[lb]{{$Q_{ai}$}}}
\put(2026,-3736){\makebox(0,0)[lb]{{$h_i$}}}
\put(2026,-3136){\makebox(0,0)[lb]{{$h_s$}}}
\put(4801,-4186){\makebox(0,0)[lb]{{$T_0$}}}
\put(4801,-2986){\makebox(0,0)[lb]{{$T_3$}}}
\put(4801,-3286){\makebox(0,0)[lb]{{$T_2$}}}
\put(4801,-3736){\makebox(0,0)[lb]{{$T_1$}}}
\end{picture}
}
\caption{Diagram of internal ice temperatures and fluxes. The hashed
layer is the snow.}
\label{fflux}
\end{figure}
Figure \ref{fflux} shows the locations of the ice and snow temperatures
and the heat fluxes. The temperature profile is assumed to be linear
between adjacent temperature points. The interior of the ice contains
``brine pockets'', leading to a prognostic equation for the temperature
$T_1$.
The surface flux to the air is:
\begin{equation}
Q_{ai} = - H\!\downarrow - LE\!\downarrow -
\epsilon_s LW\!\!\downarrow -
(1 - \alpha_s) SW\!\!\downarrow + \epsilon_s \sigma (T_3+273)^4
\end{equation}
The incoming shortwave and longwave radiations are assumed to
come from an atmospheric model. The formulae for sensible heat,
latent heat, and outgoing longwave radiation are the same as in
\citet{Parkinson} and are shown in Appendix
\ref{shortwave}. The sensible heat is a function of $T_3$, as is the
heat flux through the snow $Q_s$. Setting $Q_{ai} = Q_s$, we can solve
for $T_3$ by setting $T_3^{n+1} = T_3^n + \Delta T_3$ and linearizing in
$\Delta T_3$. As in Parkinson and Washington, if $T_3$
is found to be above the melting temperature, it is set to $T_{\rm melt}$
and the extra energy goes into melting the snow or ice:
\begin{align}
W_{ai} & = \frac{Q_{ai} - Q_{i2} }{ \rho_o L_3} \label{eqWai} \\
L_3 & \equiv \left[ E(T_3,1) - E(T_1, r) \right]
\end{align}
Note that $L_3 = (1-r)L_i$ plus a small sensible heat correction.
If there is no snow, there is an option to allow some of the incoming
shortwave radiation ($I_o = 0.17(1-\alpha_i)SW\!\!\downarrow$) to
enter the ice and contribute to heating up the ice rather than to
melting the surface ice \citep{Maykut71}. It essentially reduces $Q_{i2}$
by the amount $I_o$ in both equations (\ref{eqWai}) and (\ref{eqdEdt}).
%We can store water on the surface in melt pools to a fixed
%depth---everything extra melted at the surface is assumed to flow into
%the ocean as $W_{ro}$. The melt ponds however do not contribute to
%the heat flux computation.
Inside the ice there are brine pockets in which there is salt water
at the {\it in situ} freezing temperature. It is assumed that the ice
has a uniform overall salinity of $S_i$ and that the freezing
temperature is a linear function of salinity. The brine fraction $r$ is
given by
$$
r = \frac{S_i m }{ T_1}
$$
The enthalpy of the combined ice/brine system is given by
\begin{equation}
E(T,r) = r(L_i + C_{po}T) + (1-r) C_{pi} T
\end{equation}
Substituting in for $r$ and differentiating gives:
\begin{equation}
\frac{\partial E }{ \partial T_1} = - \frac{S_i m L_i }{ T_1^2} + C_{pi}
\end{equation}
Inside the snow, we have
\begin{equation}
Q_s = \frac{k_s }{ h_s} (T_2 - T_3)
\end{equation}
The heat conduction in the upper part of the ice layer is
\begin{equation}
Q_{I2} = \frac{ 2 k_i }{ h_i} (T_1 - T_2)
\label{qi2}
\end{equation}
These can be set equal to each other to solve for $T_2$
\begin{equation}
T_2 = \frac{T_3 + C_k T_1 }{ 1 + C_k}
\end{equation}
where
$$
C_k \equiv \frac{2 k_i h_s }{ h_i k_s}.
$$
Substituting into (\ref{qi2}), we get:
\begin{equation}
Q_s = Q_{I2} = \frac{2k_i }{ h_i} \frac{(T_1 - T_3) }{ (1 + C_k)}
\label{qsnow}
\end{equation}
Note that in the absence of snow, $C_k$ becomes zero and we recover the
formula for the no-snow case in which $T_3 = T_2$.
At the bottom of the ice, we have
\begin{equation}
Q_{I0} = \frac{2 k_i }{ h_i} (T_0 - T_1)
\end{equation}
The difference between $Q_{I0}$ and $Q_{I2}$ goes into the enthalpy of
the ice:
\begin{equation}
\rho_i h_i \left[ \frac{\partial E }{ \partial t} + \vec{v} \cdot
\nabla E \right] = Q_{I0} - Q_{I2}
\label{eqdEdt}
\end{equation}
We can use the chain rule to obtain an equation for timestepping $T_1$:
\begin{equation}
\rho_i h_i \frac{\partial E }{ \partial T_1}
\left[ \frac{\partial T_1 }{ \partial t} + \vec{v} \cdot
\nabla T_1 \right] = Q_{I0} - Q_{I2}
\end{equation}
where
\begin{align*}
Q_{I0} - Q_{I2} & = \frac{2 k_i }{ h_i} \left[ (T_0 - T_1) -
\frac{(T_1 - T_3) }{ 1 + C_k} \right] \\
& = \frac{2 k_i }{ h_i} \left[ (T_0 +
\frac{T_3 - (2 + C_k) T_1 }{ 1 + C_k} \right]
\end{align*}
When the snow gets thick enough, it pushes the ice down to the point
where seawater floods onto the ice where it can possibly refreeze.
The condition for such thick snow is $h^\star > 0$, where
\begin{equation}
h^\star = h_s - \frac{(\rho_w - \rho_i)}{\rho_s} h_i
\end{equation}
An optional model process is to convert the thick snow into ice
directly:
\begin{align}
h_s & = h_s - \frac{\rho_i h^\star}{\rho_w} \\
h_i & = h_i + \frac{\rho_s h^\star}{\rho_w}
\end{align}
\subsubsection{Ocean surface boundary conditions}
The ocean receives surface stresses from both the atmosphere and the
ice, according to the ice concentration:
\begin{align}
K_m \frac{\partial u_w }{ \partial z} & = \frac{A }{ \rho_o} \tau_{io}^x
+ \frac{1-A }{ \rho_o} \tau_{ao}^x \\
K_m \frac{\partial v_w }{ \partial z} & = \frac{A }{ \rho_o} \tau_{io}^y
+ \frac{1-A }{ \rho_o} \tau_{ao}^y
\end{align}
where the relevant variables are in Table \ref{tobc}.
\begin{table}
\hspace{35 mm}
\vbox{
\begin{tabular}{|c|c|l|} \hline
Variable & Value & Definition \\ \hline
$b$ & 3.14 & factor \\
\.E && evaporation \\
$k$ & 0.4 & von Karman's constant \\
$K_m$ && vertical viscosity of seawater \\
$\nu$ & $1.8 \times 10^{-6} m^2 s^{-1}$ &
kinematic viscosity of seawater \\
\.P && precipitation \\
$Pr$ & 13.0 & molecular Prandtl number \\
$P_{rt}$ & 0.85 & turbulent Prandtl number \\
$S$ && internal ocean salinity \\
$S_0$ && surface salinity \\
$Sc$ & 2432.0 & molecular Schmidt number \\
$\tau_{io}$ && stress on the ocean from the ice \\
$\tau_{ao}$ && stress on the ocean from the wind \\
$T$ && internal ocean temperature \\
$T_0$ && surface temperature \\
$u_\tau$ && friction velocity $|\tau_{io}|^{1/2} \rho_o^{-1/2}$ \\
$z_0$ && roughness parameter \\
\hline
\end{tabular}
}
\caption{Ocean surface variables}
\label{tobc}
\end{table}
The surface ocean is assumed to be at the freezing temperature for the
surface salinity ($T_0 = mS_0$) in the presense of ice. We also have $T$ and $S$ at the
uppermost computed ocean point ${1 \over 2} dz$ below the surface.
In order to solve for $T_0$ and $S_0$, we assume a \citep{Yaglom_1974}
logarithmic boundary layer. The upper ocean heat flux is:
\begin{equation}
\frac{F_T }{ \rho_o C_{po}} = -C_{T_z} (T_0 - T) \qquad z \rightarrow 0
\end{equation}
where
\begin{gather}
C_{T_z} = \frac{u_\tau }{ P_{rt} k^{-1}\ln (-z/z_0) + B_T} \\
B_T = b \left(\frac{z_0 u_\tau }{ \nu} \right) ^{1/2} Pr^{2/3}
\end{gather}
Likewise, we have the following equation for the surface salt flux:
\begin{equation}
F_S = -C_{S_z} (S_0 - S) \qquad z \rightarrow 0
\end{equation}
where
\begin{gather}
C_{S_z} = \frac{u_\tau }{ P_{rt} k^{-1}\ln (-z/z_0) + B_S} \\
B_S = b \left(\frac{z_0 u_\tau }{ \nu} \right) ^{1/2} Sc^{2/3}
\end{gather}
The ocean model receives the following heat and salt fluxes:
\begin{align}
F_T & = A Q_{io} + (1 - A) Q_{ao} - W_o L_o \\
F_S & = (W_o - A W_{ro}) (S_i-S_0) + (1-A)S_o (\mbox{\.P}-\mbox{\.E}) \\
W_o & \equiv A W_{io} + (1-A) W_{ao}
\end{align}
\citep{Mellor89} describe solving simultaneously for
the five unknowns $W_o$, $T_0$, $S_0$, $F_T$ and $F_S$. Instead, we use
the old value of $T_0$ to find $W_{io}$ and therefore $W_o$. Using the
new value of $W_{io}$, solve for a new value of $S_0$ and then find the
new $T_0$ as the freezing temperature for that salinity:
\begin{align}
W_{io} &= {1 \over L_o} \left[ {Q_{io} \over \rho_o} + C_{po}
C_{T_z} (T_o - T) \right] \\
S_0 &= { C_{S_z} S + (W_{ro} - W_{io})) S_i \over
C_{S_z} + W_{ro} - W_{io} }
\end{align}
\subsubsection{Frazil ice formation}
\label{frazil}
Following \citet{Steele89}, we check to see if any
of the ocean temperatures are below freezing at the end of each
timestep. If so, frazil ice is formed, changing the local
temperature and salinity. The ice that forms is assumed to
instantly float up to the surface and add to the ice layer there.
We balance the mass, heat, and salt before and after the ice
is formed:
\begin{align}
m_{w_1} & = m_{w_2} + m_i \\
m_{w_1} ( C_{pw} T_1 + L) & =
m_{w_2} (C_{pw} T_2 + L) + m_i C_{pi} T_2 \\
m_{w_1} S_1 & = m_{w_2} S_2 .
\end{align}
The variables are defined in Table \ref{frazvar}.
\begin{table}[t]
\hspace{35 mm}
\vbox{
\begin{tabular}{|c|c|l|} \hline
Variable & Value & Definition \\ \hline
$C_{pi}$ & 1994 J kg$^{-1}$ K$^{-1}$ & specific heat of ice \\
$C_{pw}$ & 3987 J kg$^{-1}$ K$^{-1}$ & specific heat of water \\
$\gamma$ & $m_i/m_{w_2}$ & fraction of water that froze \\
$L$ & 3.16e5 J kg$^{-1}$ & latent heat of fusion \\
$m_i$ && mass of ice formed \\
$m_{w_1}$ && mass of water before freezing \\
$m_{w_2}$ && mass of water after freezing \\
$m$ & $-0.0543$ & constant in freezing equation \\
% $n$ & $7.59 \times 10^{-4}$ & constant in freezing equation \\
$S_1$ && salinity before freezing \\
$S_2$ && salinity after freezing \\
$T_1$ && temperature before freezing \\
$T_2$ && temperature after freezing \\
\hline
\end{tabular}
}
\caption{Frazil ice variables}
\label{frazvar}
\end{table}
Defining $\gamma = m_i / m_{w_2}$ and dropping terms of order $\gamma^2$
leads to:
\begin{align}
T_2 & = T_1 + \gamma \left[ \frac{L }{ C_{pw}} + T_1 \left( 1
- \frac{C_{pi} }{ C_{pw}} \right) \right] \label{t2eq} \\
S_2 & = S_1 (1 + \gamma) \label{s2eq} .
\end{align}
We also want the final temperature and salinity to be on the freezing
line, which we approximate as:
\begin{equation}
T_f = m S .
% T_f = m S + n z .
\label{eq_t_f}
\end{equation}
We can then solve for $\gamma$:
\begin{equation}
% \gamma = \frac{-T_1 + mS_1 + nz }{ {L }{ C_{pw}}+ T_1 \left( 1
\gamma = \frac{-T_1 + mS_1 }{ {L }{ C_{pw}}+ T_1 \left( 1
- \frac{C_{pi} }{ C_{pw}} \right) - mS_1} .
\label{eq_gam}
\end{equation}
The ocean is checked at each depth $k$ and at each timestep for
supercooling. If the water is below freezing, the temperature and
salinity are adjusted as in equations (\ref{t2eq}) and (\ref{s2eq})
and the ice above is thickened by the amount:
\begin{equation}
\Delta h = \gamma_k \Delta z_k \frac{\rho_w}{\rho_i} .
\end{equation}
Note that \citet{Steele89} include a compressibility term in equations
(\ref{eq_t_f}) and (\ref{eq_gam}), but we use potential temperature
instead of {\em in situ} temperature.
\subsubsection{Differences from Mellor and Kantha}
We have tried to modify the \code{hakkis} model to more closely follow
\citet{Mellor89}. However, there are also ways in which we have deviated from it.
\begin{itemize}
\item Add advection of snow.
\item Add lateral melting of snow when ice is melting laterally.
% \item We iterate on the solution of $T_3$.
% \item We took a shortcut in the solution of $S_0, T_0$ for the
% surface heat and salt fluxes. We also apply them differently to the
% ocean model.
\item Add various limiters:
\begin{itemize}
\item Ice concentration:$A_{\min} \leq A \leq 1.0$, $A_{\min} =
1.e-30$.
\item Ice thickness: $h_i \geq 0.0$.
\item Snow thickness: $h_s \geq 0.0$.
\item Brine fraction: $r \leq r_{\max}$, $r_{\max} = 0.2$
\end{itemize}
\end{itemize}