- Email: [email protected]
Numerical simulations of sea ice with different advection schemes
372
2011,23(3):372-378 DOI: 10.1016/S1001-6058(10)60125-4
NUMERICAL SIMULATIONS OF SEA ICE WITH DIFFERENT ADVECTION SCHEMES* LIU Xi-ying State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China Institute of Meteorology, PLA University of Science and Technology, Nanjing 211101, China, E-mail: [email protected] (Received November 16, 2010, Revised April 4, 2011) Abstract: Numerical simulations are carried out for sea ice with four different advection schemes to study their effects on the simulation results. The sea ice model employed here is the Sea Ice Simulator (SIS) of the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model version 4b (MOM4b) and the four advection schemes are, the upwind scheme originally used in the SIS, the Multi-Dimensional Positive Advection (MDPA) scheme, the Incremental Remapping Scheme (IRS) and the Two Step Shape Preserving (TSSP) scheme. The latter three schemes are newly introduced. To consider the interactions between sea ice and ocean, a mixed layer ocean model is introduced and coupled to the SIS. The coupled model uses a tri-polar coordinate with 120×65 grids, covering the whole earth globe, in the horizontal plane. Simulation results in the northern high latitudes are analyzed. In all simulations, the model reproduces the seasonal variation of sea ice in the northern high latitudes well. Compared with the results from the observation, the sea ice model produces some extra sea ice coverage in the Greenland Sea and Barents Sea in winter due to the exclusion of ocean current effects and the smaller simulated sea ice thickness in the Arctic basin. There are similar features among the results obtained with the introduced three advection schemes. The simulated sea ice thickness with the three newly introduced schemes are all smaller than that of the upwind scheme and the simulated sea ice velocities of movement are all smaller than that of the upwind scheme. There are more similarities shared in the results obtained with the MPDA and TSSP schemes. Key words: advection scheme, sea ice, simulation
Introduction As an important component of cryosphere, the sea ice plays an important role in the global climate system through ocean-sea ice-atmosphere interactions. In the simulation of a climate system, an accurate representation of the sea ice has an important impact on the overall results. It was found that the simulation results of climate in high latitudes are strongly influenced by the disposition of the sea ice[1,2]. Particularly, in the simulation with an ocean-atmosphere coupled model, the behavior of the sea ice seems * Project supported by the National Natural Science Foundation of China (Grant No. 40876101), the National High Technology Research and Development Program of China (863 Program, Grant No. 2010AA012304). Biography: LIU Xi-ying (1970-), Male, Ph. D., Associate Professor
to be critical to the results[3]. So, it is a challenge to represent the sea ice appropriately for the development of climate models. It is also shown from numerical experiments that the processes related with the sea ice dynamics play an important role in climate. For example, in a long period simulation experiment, the wind-driven changes in the sea-ice area are about twice as large as those due to thermodynamic (i.e., radiative) forcing[4]. Hence, the sensibility of a sea ice model to the variation of the thermal forcing is larger if only the thermodynamic forcing is considered than if both thermodynamic forcing and sea ice dynamics are considered. Similarly, when the dynamical variation of the sea ice is included in the coupled climate system modeling, the sensitivity of an atmospheric general circulation model to the global warming effect due to the doubling of CO2 is decreased[5]. Based on the results from a comparison of the simulated sea ice volume with the observed data by submarine in the
373
northern Arctic, it is concluded that to reveal the interannual variability and to better understand the real climatic signals, an appropriate sea ice dynamics model is needed[6]. In Arctic climate simulations, the sea ice is found to vary on a large scale. On a small-meso scale, the sea ice dynamic characteristics are quite different from those on a large scale[7]. When the sea ice dynamics is considered in the simulation, the ice velocity as well as the corresponding advection of the ice masses can be calculated. The accurate and efficient representation of advection is a key element in any numerical model for the fluid motion[8,9], which is true for the sea ice simulation as well. Under the assumption of continuum, the algorithm used to solve the advection process numerically for the sea ice variations is closely related to the quality of the solutions and the numerical scheme for advection has a profound impact on the simulated sea ice advection and hence on the overall results. In the sea ice simulation, the variations of the sea ice area ( A ) and the sea ice volume ( V ) can be described by[10]
the sea ice internal interaction, and W a and W o are the stresses related with atmosphere and ocean, respectively. In the numerical simulation of the sea ice, the horizontal transport of the sea ice is usually represented by the item <(uFn ) (see Eq.(1) and Eq.(2)). Here Fn can be An , Vn and others. It is obvious that the chosen computation scheme may have an importantl influence on the simulation results. In this study, three new advection schemes are introduced into the sea ice model and numerical simulations of the sea ice with different advection schemes are performed to study their effects on simulation results. The four advection schemes, including the three newly introduced ones and the original one in the model, will be discussed in Section 1. The experiment configurations as well as simulation results will be discussed in Section 2 and conclusions will be given in Section 3.
where ST is the source/sink due to thermodynamical
1. Sea ice advection schemes 1.1 Upwind scheme The sea ice model Sea Ice Simulator (SIS) in Modular Ocean Model, version 4 (MOM4) of Geophysical Fluid Dynamics Laboratory (GFDL)[11,12] deals with the sea ice advection with the upwind scheme. In the upwind scheme, the one-dimensional advection equation
processes and the transport in thickness and S M is the source/sink due to the mechanical redistribution. The thickness of the sea ice can be calculated by hn =
w\ w = u\ wt wx
wAn = STAn < uAn + S MAn wt
(1)
wVn = STVn < uVn + S MVn wt
(2)
Vn / An . The subscript n represents the category of the sea ice thickness and u is the velocity vector of the sea ice current. Similar equations for the internal energy of the sea ice, the snow volume and the sea ice surface temperature weighted by the sea ice area can be established. In the equation of momentum describing the movement of the sea ice, since the sea ice velocity is small, the item of the nonlinear advection of momentum is small compared to other items and is usually neglected[10] and the equation takes the form m
wu = mfk u u + W a + W o + mg H o +
(3)
where m is the mass of the sea ice and the snow (if it is on the sea ice surface), f is the Coriolis parameter, k is the unit vector in the vertical direction, g is the gravity acceleration, H o is the height of the sea surface, V is the stress tensor related with
(4)
is discretized as
\ im +1 = \ im ª¬ F \ im ,\ im+1 , uim+1/ 2 F \ im1 ,\ im , uim1/ 2 º¼ (5) where
F \ i ,\ i +1 , u = ª¬ u + u \ i + u u \ i +1 º¼
't 2'x
The disadvantage of the upwind scheme lies in its strong implicit diffusion. 1.2 Introduction of the multi-dimensional positive advection scheme To diminish the strong diffusion in the upwind scheme, a Multi-Dimensional Positive Advection (MDPA) scheme is introduced into the SIS. The basic idea of the MDPA is to use an imported antidiffusion velocity to correct the result of the upwind scheme. The scheme is as follows[13] (1) Get an initial guess value \ with the up-
374
wind scheme
If the relation (9) does not hold true, the upwind scheme is used, otherwise \ i is corrected with formula (10).
\ i = \ im ª¬ F \ im ,\ im+1 , uim+1/ 2 F \ im1 ,\ im , uim1/ 2 º¼ (6) (2) Correct the \ diffusion velocity
\
with an introduced anti-
= \ ª¬ F \ ,\ , u c
i
m +1 i
i
i +1
m i 1/ 2
F \
* i 1
,\ , u c * i
m i 1/ 2
't ªuim+1/ 2 \ im+1 + \ im uim1/ 2 < 2'x ¬
\ im +1 = \ im
\
º¼ (7)
u
m i
+ \ im1 º¼ +
m i 1/ 2
\
m i
The anti-diffusion velocity is defined as m uic+1/ 2
u2
m i
n
't ªuim+1/ 2 \ im+1 + \ im uim1/ 2 < 2'x ¬
+ \ im1 ¼º +
m i 1/ 2
't ª 2 m 't m u \ i +1 \ im 2'x ¬« i +1/ 2 'x
't m º \ i \ im1 » 'x ¼
(8)
(2) Correct Determine if \ i meets the criterion
\ immin d \ i d \ immin
(9)
where
\
m i min
(10)
n
u i ±1/ 2 = Ci r1/ 2 u i r1/ 2 ,
The Sstep (2) can be iterated several times. As long as the velocity satisfies the stability condition, the subsequent anti-diffusion velocity in the iteration would be diminished gradually. Following Smolarkiewicz’s advice[13], we use Eq.(7) three times in the simulation experiments. 1.3 Introduction of the two step shape preserving scheme A Two Step Shape Preserving (TSSP) scheme, which performs well in the Arctic Ocean simulation[14] is also introduced into the SIS. The TSSP scheme also involves two steps, guess with the Lax-Wendroff scheme and correct with the upwind scheme. As their combination, the strong points of the two schemes are kept. (1) Guess with the Lax-Wendroff scheme
\
\ im1 º ¼
where
0.5 uim+1/ 2 'x uim+1/ 2 2 't w\ = \ wx
\ i = \ im
't ª m u \ im+1 \ im 2'x ¬ i +1/ 2
Ci r1/ 2 = Ci r1/ 2 + 1 Ci r1/ 2
n
u i r1/ 2 't
'x
,
§ A + Ai Ai r1 + Ai r1 · Ci r1/ 2 = 0.5 ¨ i + ¸ ¨ A +H Ai r1 + H ¸¹ © i 0.25
A
i
+ Ai Ai r1 + Ai r1 Ai Ai r1 + H
Here Ai = (\ i \ inmax ) (\ i \ inmin ) , and H is a small
quantity used to avoid the case of dividing by zero. 1.4 Introduction of the incremental remapping scheme The Incremental Remapping Scheme (IRS), proposed by Dukowicz and Baumgardner[15] and used in the sea ice simulation for the first time by Lipscomb and Hunke[16], is introduced into the SIS. The algorithm of the IRS can be summarized as follows: (1) Reconstruct sea ice area and tracer fields Given mean values of the ice area and tracer fields in each grid cell, construct linear approximations of these fields and limit the field gradients to preserve monotonicity. Take the advection of the sea ice area as an example, the sea ice area (denoted by a ) in the center of the position (i, j ) (denoted by subscript c ) can be approximated linearly by
ac = a ax x a y y where
= min \ ,\ ,\ m i 1
m i
m i +1
,
\ immin = max \ im1 ,\ im ,\ im+1
a=³
A
adA wa wa , ay = D , , ax = D wx wy A
375
x=³
A
ydA xdA , y=³ . A A A
D = min(D max , D min )
is used to limit the field
gradients to preserve monotonicity. Here D max = (amax a ) /(amax a ) , D min = (amin a ) / (amin a ) ,
amax and amin are the maximum and the minimum of a at point (i, j ) and nearby place. amax and amin are the maximum and the minimum of a at point (i, j ) and nearby place after reconstruction. If amax ! amax or amin amin , the monotonicity is destroyed and the coefficient D is used. (2) Locate departure triangles Given ice velocities at grid cell corners, identify departure regions for the fluxes across each cell edge. Divide these departure regions into triangles and compute the coordinates of the triangle vertices. (3) Integrate fluxes Integrate these fields over the departure triangles to obtain the area, volume, and the energy fluxes across each cell edge. (4) Update state variables Transfer the fluxes across cell edges and update the state variables. For details of calculation, the reader can refer to [16]. 2. Numerical simulations with four different advection schemes 2.1 Description of simulation experiments Four numerical experiments with different advection schemes are designed. These experiments are identical except for the advection schemes. A mixed layer ocean model, which is used in Community Sea Ice Model, version 4 (CSIM4) is introduced into SIS to calculate variations of the Sea Surface Temperature (SST). The depth of the mixed layer in the ocean model is 30 m. The salinity and the initial SST are taken from the dataset Polar Science Center Hydrographic Climatology (PHC)[17]. The atmospheric forcing data come from the dataset ERA40[18] and 12 monthly means averaged over 44 years (1958-2001) are used. Except for the radiation fluxes, the atmospheric state quantities are supplied to the SIS and the fluxes of momentum, the sensible heat and the latent heat are calculated in the model. In the experiments, the effects of the sea surface current and the tilt of the sea surface are not considered. The effect of the precipitation is also not considered. The model uses a tripolar coordinate with 120×65 grids in the horizontal plane. 2.2 Results The model’s results are integrated for 45 years for all configurations. The ice volume grows inten-
sively in the first 5 years. Then its growth rate diminishes and the total ice volume tends to stabilize in a few years (the figure not shown). To exclude the influence of the adjustment processes in the early integrations, the results of the last 15 years are to be analyzed. Since the sea ice extent reaches its maximum in February-March and its minimum in AugustSeptember in the northern hemisphere, March and September will be used to represent the winter and the summer, respectively. 2.2.1 Results with the upwind scheme In winter, the sea ice distribution in the Hudson Bay, Fox Bay, Baffin Bay, Bering Sea, Okhotsk Sea, Davies strait and Laborador Sea are all reproduced (see Fig.1(a)). Compared with the results of observation[19], there are some extra sea ice cover in Greenland Sea and Barents Sea. Influenced by North Atlantic warm current, there is no sea ice south of the Spitsbergen Island and there is no sea ice in the southern part of Barents Sea as well. These features cannot be reproduced due to the neglect of the effect of the ocean current in the model. In summer, the sea ice extent shrinks northward. Sea ices in the Hadson Bay, Bering Sea, Okhotsk Sea, Davies strait and Laborador Sea are melted away. These features are reproduced well (Fig.1(b)). But there are some extra ices in Kara Sea and North of Baffin Bay in simulation results. The simulated sea ice thickness is smaller than the observation value[20]. These discrepancies are related with atmospheric forcing as well as model parameters.
Fig.1 Simulated distribution of sea ice concentration in March (a) and September (b) with upwind scheme
The simulated sea ice current is plotted after transforming the result from the tri-polar coordinate to the rotated spherical coordinate with the geographical North Pole as the center in the model equator. It is shown that the sea ice flow from the Greenland Sea to the North Atlantic Ocean is reproduced but the model cannot reproduce the Beaufort Gyre and the Transpolar Drift well (the figure not shown). Simulation results of the sea ice current are related closely with atmospheric wind forcing. The failure in reproducing the Beaufort Gyre may also be due to the low model resolution used.
376
2.2.2 Results with the MDPA scheme Compared with the results of the upwind scheme, the area covered by the sea ice with concentration larger than 0.9 shrinks northward in winter. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of the Okhotsk Sea becomes larger (see Fig.3(a) and Fig.1(a)). In summer, the sea ice extent with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea shrinks northward (see Fig.3b and Fig.1(b)). The sea ice thickness simulated with the MDPA is smaller than that with the upwind scheme in the center of the Arctic Ocean (see Fig.4 and Fig.2). There is little difference in the sea ice current between the results of the two schemes in winter. But the velocity becomes smaller in most area in summer compared to the results with the upwind scheme (the figure not shown).
2.2.3 Results with the TSSP scheme Similar features can be found between the results of the TSSP and MDPA schemes (comparing Fig.5 with Fig.3, and Fig.6 with Fig.4). As compared with the results of the upwind scheme, the simulated sea ice area obtained with the TSSP scheme with concentration larger than 0.9 shrinks northward in winter. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of Okhotsk Sea becomes larger (see Fig.5(a) and Fig.1(a)). In summer, the sea ice cover with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea shrinks northward (see Fig.5(b) and Fig.1(b)). The sea ice thickness simulated with the TSSP scheme is smaller than that of the upwind scheme in the center of the Arctic Ocean (see Fig.6 and Fig.2). There is little difference in the sea ice current between the results obtained with the TSSP and upwind schemes in winter. But the velocity becomes smaller in most area in summer as compared to the results obtained with the upwind scheme (the figure not shown).
Fig.2 As in Fig.1 except for sea ice thickness
Fig.5 As in Fig.1 except for TSSP scheme
Fig.3 As in Fig.1 except for MDPA scheme
Fig.6 As in Fig.2 except for TSSP scheme
Fig.4 As in Fig.2 except for MDPA scheme
2.2.4 Results with the IRS scheme The simulated results obtained with the IRS scheme are shown in Fig.7 and Fig.8. Compared to the results of the upwind scheme, there is a larger area with sea ice concentration over 0.9 in Bering Sea and Okhotsk Sea whereas there is a smaller area with sea ice concentration over 0.9 in Labrador Sea in winter
377
(see Fig.7(a) and Fig.1(a)). In summer, the ice covered area with sea ice concentration over 0.9 in Chukchi Sea shrinks northward and the ice concentration is bigger in places near the continent in the Greenland Sea (see Fig.7(b) and Fig.1(b)). Compared to the results of the upwind scheme, the simulated sea ice thickness is smaller (See Fig.8 and Fig.2) and the simulated sea ice velocity is smaller in most ice covered area (the figure not shown).
velocities of movement are all smaller than that of the upwind scheme. There are more similarities shared in the results obtained with the MDPA and TSSP schemes. In both schemes, the simulated sea ice extent with concentration larger than 0.9 is smaller in winter as compared to results of the upwind scheme. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of Okhotsk Sea becomes larger. In summer, the sea ice extent with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea contracts northward. Although the performances of different advection schemes can be analyzed in simple cases, in real complicated applications, due to influences of other mechanisms and complicated interactions involved, more uncertainties may arise. References
Fig.7 As in Fig.1 except for IRS scheme
[1]
[2]
[3] [4]
Fig.8 As in Fig.2 except for IRS scheme
3. Conclusion In the study, we introduce three new advection schemes into the sea ice model SIS. With an introduction of a simple mixed layer ocean model into the SIS to calculate the variations of the sea surface temperature, four numerical simulations of the sea ice using different advection schemes are performed. Simulation results in the northern high latitudes are analyzed. It is shown that the SIS can reproduce the seasonal variations of the sea ice well with whichever advection scheme. In all simulations, the SIS produces some extra sea ice cover in the Greenland Sea and Barents Sea in winter due to the neglect of the influence of ocean currents. The sea ice model produces smaller ice thickness in the Arctic basin than the observation value, which may be due to an inappropriate parameterization of the heat flux from the lower ocean water. There are similar features among results obtained with the introduced three advection schemes. The simulated sea ice thicknesses are all smaller than that of the upwind scheme and the simulated sea ice
[5] [6] [7]
[8]
[9]
[10] [11] [12]
PARKINSON C. L., RIND D. and HEALY R. J. et al. The impact of sea ice concentration on climate model simulations with the GISS GCM[J]. J. Climate, 2001, 14(12): 2606-2623. ALEXANDER M. A., BHATT U. S. and WALSH J. E. et al. The atmospheric response to realistic Arctic sea ice anomalies in an AGCM during winter[J]. J. Climate, 2004, 17(5): 890-905. HOLLAND M., BIT C. and HUNKE E. et al. Influence of the sea ice thickness distribution on polar climate in CCSM3[J] J. Climate, 2006, 19(11): 2398-2414. SEDLÁCEK J., MYSAK L. A. Sensitivity of sea ice to wind-stress and radiative forcing since 1500: A model study of the little ice age and beyond[J]. Clim. Dyn., 2009, 32(6): 817-831. HOLLAND M. M., BITZ C. M. Polar amplification of climate change in coupled models[J]. Clim. Dyn., 2003, 21(3-4): 221-232 HILMER M., LEMKE P. On the decrease of arctic sea ice volume[J]. Geophys. Res. Let., 2000, 27(22): 37513754. WANG Gang, JI Shun-ying and LV He-xiang et al. Drucker-prager yield criteria in viscoelastic-plastic constitutive model for the study of sea ice dynamics[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(6): 714722. WANG H., SKAMAROCK W. C. and FEINGOLD G. Evaluation of scalar advection schemes in the advanced research WRF model using large-eddy simulations of aerosol-cloud interactions[J]. Mon. Wea. Rev., 2009, 137(8): 2547-2558. YANG Wei, LIU Shu-hong and WU Yu-lin. An unsplit lagrangian advection scheme for volume of fluid method[J]. Journal of Hydrodynamics, 2010, 22(1): 73-80. MARTIN T. Arctic sea ice dynamics: Drift and ridging in numerical models and observations[D]. Ph. D. Thesis, Potsdam, Germany: Alfred Wegener Institute, 2007. WINTON M. A reformulated three-layer sea ice model[J]. Journal of Atmospheric and Oceanic Technology, 2000, 17(4): 525-531. DELWORTH T. L. and Couuthors. GFDL’s CM2 global coupled climate models. Part I: Formulation and
378
[13]
[14]
[15] [16]
simulation characteristics[J]. J. Climate, 2006, 19(5): 643-674. SMOLARKIEWICZ P. K. Multidimensional positive definite advection transport algorithm: An overview[J]. International Journal for Numerical Methods in Fluids, 2006, 50(10): 1123-1144. LIU Xi-ying, ZHANG Xue-hong and YU Ru-cong et al. Fine-resolution simulation of surface current and sea ice in the Arctic Mediterranean Seas[J]. Chinese Journal of Oceanology and Limnology, 2007, 25(2): 132-138. DUKOWICZ J. K., BAUMGARDNER J. R. Incremental remapping as a transport/advection algorithm[J]. J. Comput. Phys., 2000, 160(1): 318-335. LIPSCOMB W. H., HUNKE E. C. Modeling sea ice transport using incremental remapping[J]. Mon.Wea. Rev., 2004, 132(6): 1341-1354.
[17] [18] [19]
[20]
STEELE M., MORLEY R. and ERMOLD R. PHC: A global ocean hydrography with a high quality Arctic Ocean[J]. J. Climate, 2001, 14(9): 2079-2087. UPPALA S. M., KÅLLBERG P. W. and SIMMONS A. J. et al. The ERA-40 re-analysis[J]. Quart. J. R. Meteorol. Soc., 2005, 131(612): 2961-3012. LIU X. Y., ZHANG X. H. and YU Y. Q. et al. Mean climatic characteristics in high northern latitudes in an ocean-sea ice-atmosphere coupled model[J]. Adv. Atmos. Sci., 2004, 21(2): 236-244. BOURKE R. H., GARRETT R. P. Sea-ice thickness distribution in the Arctic Ocean[J]. Cold Reg. Sci. Technol., 1987, 13(3): 259-280.
2011,23(3):372-378 DOI: 10.1016/S1001-6058(10)60125-4
NUMERICAL SIMULATIONS OF SEA ICE WITH DIFFERENT ADVECTION SCHEMES* LIU Xi-ying State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China Institute of Meteorology, PLA University of Science and Technology, Nanjing 211101, China, E-mail: [email protected] (Received November 16, 2010, Revised April 4, 2011) Abstract: Numerical simulations are carried out for sea ice with four different advection schemes to study their effects on the simulation results. The sea ice model employed here is the Sea Ice Simulator (SIS) of the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model version 4b (MOM4b) and the four advection schemes are, the upwind scheme originally used in the SIS, the Multi-Dimensional Positive Advection (MDPA) scheme, the Incremental Remapping Scheme (IRS) and the Two Step Shape Preserving (TSSP) scheme. The latter three schemes are newly introduced. To consider the interactions between sea ice and ocean, a mixed layer ocean model is introduced and coupled to the SIS. The coupled model uses a tri-polar coordinate with 120×65 grids, covering the whole earth globe, in the horizontal plane. Simulation results in the northern high latitudes are analyzed. In all simulations, the model reproduces the seasonal variation of sea ice in the northern high latitudes well. Compared with the results from the observation, the sea ice model produces some extra sea ice coverage in the Greenland Sea and Barents Sea in winter due to the exclusion of ocean current effects and the smaller simulated sea ice thickness in the Arctic basin. There are similar features among the results obtained with the introduced three advection schemes. The simulated sea ice thickness with the three newly introduced schemes are all smaller than that of the upwind scheme and the simulated sea ice velocities of movement are all smaller than that of the upwind scheme. There are more similarities shared in the results obtained with the MPDA and TSSP schemes. Key words: advection scheme, sea ice, simulation
Introduction As an important component of cryosphere, the sea ice plays an important role in the global climate system through ocean-sea ice-atmosphere interactions. In the simulation of a climate system, an accurate representation of the sea ice has an important impact on the overall results. It was found that the simulation results of climate in high latitudes are strongly influenced by the disposition of the sea ice[1,2]. Particularly, in the simulation with an ocean-atmosphere coupled model, the behavior of the sea ice seems * Project supported by the National Natural Science Foundation of China (Grant No. 40876101), the National High Technology Research and Development Program of China (863 Program, Grant No. 2010AA012304). Biography: LIU Xi-ying (1970-), Male, Ph. D., Associate Professor
to be critical to the results[3]. So, it is a challenge to represent the sea ice appropriately for the development of climate models. It is also shown from numerical experiments that the processes related with the sea ice dynamics play an important role in climate. For example, in a long period simulation experiment, the wind-driven changes in the sea-ice area are about twice as large as those due to thermodynamic (i.e., radiative) forcing[4]. Hence, the sensibility of a sea ice model to the variation of the thermal forcing is larger if only the thermodynamic forcing is considered than if both thermodynamic forcing and sea ice dynamics are considered. Similarly, when the dynamical variation of the sea ice is included in the coupled climate system modeling, the sensitivity of an atmospheric general circulation model to the global warming effect due to the doubling of CO2 is decreased[5]. Based on the results from a comparison of the simulated sea ice volume with the observed data by submarine in the
373
northern Arctic, it is concluded that to reveal the interannual variability and to better understand the real climatic signals, an appropriate sea ice dynamics model is needed[6]. In Arctic climate simulations, the sea ice is found to vary on a large scale. On a small-meso scale, the sea ice dynamic characteristics are quite different from those on a large scale[7]. When the sea ice dynamics is considered in the simulation, the ice velocity as well as the corresponding advection of the ice masses can be calculated. The accurate and efficient representation of advection is a key element in any numerical model for the fluid motion[8,9], which is true for the sea ice simulation as well. Under the assumption of continuum, the algorithm used to solve the advection process numerically for the sea ice variations is closely related to the quality of the solutions and the numerical scheme for advection has a profound impact on the simulated sea ice advection and hence on the overall results. In the sea ice simulation, the variations of the sea ice area ( A ) and the sea ice volume ( V ) can be described by[10]
the sea ice internal interaction, and W a and W o are the stresses related with atmosphere and ocean, respectively. In the numerical simulation of the sea ice, the horizontal transport of the sea ice is usually represented by the item <(uFn ) (see Eq.(1) and Eq.(2)). Here Fn can be An , Vn and others. It is obvious that the chosen computation scheme may have an importantl influence on the simulation results. In this study, three new advection schemes are introduced into the sea ice model and numerical simulations of the sea ice with different advection schemes are performed to study their effects on simulation results. The four advection schemes, including the three newly introduced ones and the original one in the model, will be discussed in Section 1. The experiment configurations as well as simulation results will be discussed in Section 2 and conclusions will be given in Section 3.
where ST is the source/sink due to thermodynamical
1. Sea ice advection schemes 1.1 Upwind scheme The sea ice model Sea Ice Simulator (SIS) in Modular Ocean Model, version 4 (MOM4) of Geophysical Fluid Dynamics Laboratory (GFDL)[11,12] deals with the sea ice advection with the upwind scheme. In the upwind scheme, the one-dimensional advection equation
processes and the transport in thickness and S M is the source/sink due to the mechanical redistribution. The thickness of the sea ice can be calculated by hn =
w\ w = u\ wt wx
wAn = STAn < uAn + S MAn wt
(1)
wVn = STVn < uVn + S MVn wt
(2)
Vn / An . The subscript n represents the category of the sea ice thickness and u is the velocity vector of the sea ice current. Similar equations for the internal energy of the sea ice, the snow volume and the sea ice surface temperature weighted by the sea ice area can be established. In the equation of momentum describing the movement of the sea ice, since the sea ice velocity is small, the item of the nonlinear advection of momentum is small compared to other items and is usually neglected[10] and the equation takes the form m
wu = mfk u u + W a + W o + mg H o +
(3)
where m is the mass of the sea ice and the snow (if it is on the sea ice surface), f is the Coriolis parameter, k is the unit vector in the vertical direction, g is the gravity acceleration, H o is the height of the sea surface, V is the stress tensor related with
(4)
is discretized as
\ im +1 = \ im ª¬ F \ im ,\ im+1 , uim+1/ 2 F \ im1 ,\ im , uim1/ 2 º¼ (5) where
F \ i ,\ i +1 , u = ª¬ u + u \ i + u u \ i +1 º¼
't 2'x
The disadvantage of the upwind scheme lies in its strong implicit diffusion. 1.2 Introduction of the multi-dimensional positive advection scheme To diminish the strong diffusion in the upwind scheme, a Multi-Dimensional Positive Advection (MDPA) scheme is introduced into the SIS. The basic idea of the MDPA is to use an imported antidiffusion velocity to correct the result of the upwind scheme. The scheme is as follows[13] (1) Get an initial guess value \ with the up-
374
wind scheme
If the relation (9) does not hold true, the upwind scheme is used, otherwise \ i is corrected with formula (10).
\ i = \ im ª¬ F \ im ,\ im+1 , uim+1/ 2 F \ im1 ,\ im , uim1/ 2 º¼ (6) (2) Correct the \ diffusion velocity
\
with an introduced anti-
= \ ª¬ F \ ,\ , u c
i
m +1 i
i
i +1
m i 1/ 2
F \
* i 1
,\ , u c * i
m i 1/ 2
't ªuim+1/ 2 \ im+1 + \ im uim1/ 2 < 2'x ¬
\ im +1 = \ im
\
º¼ (7)
u
m i
+ \ im1 º¼ +
m i 1/ 2
\
m i
The anti-diffusion velocity is defined as m uic+1/ 2
u2
m i
n
't ªuim+1/ 2 \ im+1 + \ im uim1/ 2 < 2'x ¬
+ \ im1 ¼º +
m i 1/ 2
't ª 2 m 't m u \ i +1 \ im 2'x ¬« i +1/ 2 'x
't m º \ i \ im1 » 'x ¼
(8)
(2) Correct Determine if \ i meets the criterion
\ immin d \ i d \ immin
(9)
where
\
m i min
(10)
n
u i ±1/ 2 = Ci r1/ 2 u i r1/ 2 ,
The Sstep (2) can be iterated several times. As long as the velocity satisfies the stability condition, the subsequent anti-diffusion velocity in the iteration would be diminished gradually. Following Smolarkiewicz’s advice[13], we use Eq.(7) three times in the simulation experiments. 1.3 Introduction of the two step shape preserving scheme A Two Step Shape Preserving (TSSP) scheme, which performs well in the Arctic Ocean simulation[14] is also introduced into the SIS. The TSSP scheme also involves two steps, guess with the Lax-Wendroff scheme and correct with the upwind scheme. As their combination, the strong points of the two schemes are kept. (1) Guess with the Lax-Wendroff scheme
\
\ im1 º ¼
where
0.5 uim+1/ 2 'x uim+1/ 2 2 't w\ = \ wx
\ i = \ im
't ª m u \ im+1 \ im 2'x ¬ i +1/ 2
Ci r1/ 2 = Ci r1/ 2 + 1 Ci r1/ 2
n
u i r1/ 2 't
'x
,
§ A + Ai Ai r1 + Ai r1 · Ci r1/ 2 = 0.5 ¨ i + ¸ ¨ A +H Ai r1 + H ¸¹ © i 0.25
A
i
+ Ai Ai r1 + Ai r1 Ai Ai r1 + H
Here Ai = (\ i \ inmax ) (\ i \ inmin ) , and H is a small
quantity used to avoid the case of dividing by zero. 1.4 Introduction of the incremental remapping scheme The Incremental Remapping Scheme (IRS), proposed by Dukowicz and Baumgardner[15] and used in the sea ice simulation for the first time by Lipscomb and Hunke[16], is introduced into the SIS. The algorithm of the IRS can be summarized as follows: (1) Reconstruct sea ice area and tracer fields Given mean values of the ice area and tracer fields in each grid cell, construct linear approximations of these fields and limit the field gradients to preserve monotonicity. Take the advection of the sea ice area as an example, the sea ice area (denoted by a ) in the center of the position (i, j ) (denoted by subscript c ) can be approximated linearly by
ac = a ax x a y y where
= min \ ,\ ,\ m i 1
m i
m i +1
,
\ immin = max \ im1 ,\ im ,\ im+1
a=³
A
adA wa wa , ay = D , , ax = D wx wy A
375
x=³
A
ydA xdA , y=³ . A A A
D = min(D max , D min )
is used to limit the field
gradients to preserve monotonicity. Here D max = (amax a ) /(amax a ) , D min = (amin a ) / (amin a ) ,
amax and amin are the maximum and the minimum of a at point (i, j ) and nearby place. amax and amin are the maximum and the minimum of a at point (i, j ) and nearby place after reconstruction. If amax ! amax or amin amin , the monotonicity is destroyed and the coefficient D is used. (2) Locate departure triangles Given ice velocities at grid cell corners, identify departure regions for the fluxes across each cell edge. Divide these departure regions into triangles and compute the coordinates of the triangle vertices. (3) Integrate fluxes Integrate these fields over the departure triangles to obtain the area, volume, and the energy fluxes across each cell edge. (4) Update state variables Transfer the fluxes across cell edges and update the state variables. For details of calculation, the reader can refer to [16]. 2. Numerical simulations with four different advection schemes 2.1 Description of simulation experiments Four numerical experiments with different advection schemes are designed. These experiments are identical except for the advection schemes. A mixed layer ocean model, which is used in Community Sea Ice Model, version 4 (CSIM4) is introduced into SIS to calculate variations of the Sea Surface Temperature (SST). The depth of the mixed layer in the ocean model is 30 m. The salinity and the initial SST are taken from the dataset Polar Science Center Hydrographic Climatology (PHC)[17]. The atmospheric forcing data come from the dataset ERA40[18] and 12 monthly means averaged over 44 years (1958-2001) are used. Except for the radiation fluxes, the atmospheric state quantities are supplied to the SIS and the fluxes of momentum, the sensible heat and the latent heat are calculated in the model. In the experiments, the effects of the sea surface current and the tilt of the sea surface are not considered. The effect of the precipitation is also not considered. The model uses a tripolar coordinate with 120×65 grids in the horizontal plane. 2.2 Results The model’s results are integrated for 45 years for all configurations. The ice volume grows inten-
sively in the first 5 years. Then its growth rate diminishes and the total ice volume tends to stabilize in a few years (the figure not shown). To exclude the influence of the adjustment processes in the early integrations, the results of the last 15 years are to be analyzed. Since the sea ice extent reaches its maximum in February-March and its minimum in AugustSeptember in the northern hemisphere, March and September will be used to represent the winter and the summer, respectively. 2.2.1 Results with the upwind scheme In winter, the sea ice distribution in the Hudson Bay, Fox Bay, Baffin Bay, Bering Sea, Okhotsk Sea, Davies strait and Laborador Sea are all reproduced (see Fig.1(a)). Compared with the results of observation[19], there are some extra sea ice cover in Greenland Sea and Barents Sea. Influenced by North Atlantic warm current, there is no sea ice south of the Spitsbergen Island and there is no sea ice in the southern part of Barents Sea as well. These features cannot be reproduced due to the neglect of the effect of the ocean current in the model. In summer, the sea ice extent shrinks northward. Sea ices in the Hadson Bay, Bering Sea, Okhotsk Sea, Davies strait and Laborador Sea are melted away. These features are reproduced well (Fig.1(b)). But there are some extra ices in Kara Sea and North of Baffin Bay in simulation results. The simulated sea ice thickness is smaller than the observation value[20]. These discrepancies are related with atmospheric forcing as well as model parameters.
Fig.1 Simulated distribution of sea ice concentration in March (a) and September (b) with upwind scheme
The simulated sea ice current is plotted after transforming the result from the tri-polar coordinate to the rotated spherical coordinate with the geographical North Pole as the center in the model equator. It is shown that the sea ice flow from the Greenland Sea to the North Atlantic Ocean is reproduced but the model cannot reproduce the Beaufort Gyre and the Transpolar Drift well (the figure not shown). Simulation results of the sea ice current are related closely with atmospheric wind forcing. The failure in reproducing the Beaufort Gyre may also be due to the low model resolution used.
376
2.2.2 Results with the MDPA scheme Compared with the results of the upwind scheme, the area covered by the sea ice with concentration larger than 0.9 shrinks northward in winter. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of the Okhotsk Sea becomes larger (see Fig.3(a) and Fig.1(a)). In summer, the sea ice extent with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea shrinks northward (see Fig.3b and Fig.1(b)). The sea ice thickness simulated with the MDPA is smaller than that with the upwind scheme in the center of the Arctic Ocean (see Fig.4 and Fig.2). There is little difference in the sea ice current between the results of the two schemes in winter. But the velocity becomes smaller in most area in summer compared to the results with the upwind scheme (the figure not shown).
2.2.3 Results with the TSSP scheme Similar features can be found between the results of the TSSP and MDPA schemes (comparing Fig.5 with Fig.3, and Fig.6 with Fig.4). As compared with the results of the upwind scheme, the simulated sea ice area obtained with the TSSP scheme with concentration larger than 0.9 shrinks northward in winter. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of Okhotsk Sea becomes larger (see Fig.5(a) and Fig.1(a)). In summer, the sea ice cover with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea shrinks northward (see Fig.5(b) and Fig.1(b)). The sea ice thickness simulated with the TSSP scheme is smaller than that of the upwind scheme in the center of the Arctic Ocean (see Fig.6 and Fig.2). There is little difference in the sea ice current between the results obtained with the TSSP and upwind schemes in winter. But the velocity becomes smaller in most area in summer as compared to the results obtained with the upwind scheme (the figure not shown).
Fig.2 As in Fig.1 except for sea ice thickness
Fig.5 As in Fig.1 except for TSSP scheme
Fig.3 As in Fig.1 except for MDPA scheme
Fig.6 As in Fig.2 except for TSSP scheme
Fig.4 As in Fig.2 except for MDPA scheme
2.2.4 Results with the IRS scheme The simulated results obtained with the IRS scheme are shown in Fig.7 and Fig.8. Compared to the results of the upwind scheme, there is a larger area with sea ice concentration over 0.9 in Bering Sea and Okhotsk Sea whereas there is a smaller area with sea ice concentration over 0.9 in Labrador Sea in winter
377
(see Fig.7(a) and Fig.1(a)). In summer, the ice covered area with sea ice concentration over 0.9 in Chukchi Sea shrinks northward and the ice concentration is bigger in places near the continent in the Greenland Sea (see Fig.7(b) and Fig.1(b)). Compared to the results of the upwind scheme, the simulated sea ice thickness is smaller (See Fig.8 and Fig.2) and the simulated sea ice velocity is smaller in most ice covered area (the figure not shown).
velocities of movement are all smaller than that of the upwind scheme. There are more similarities shared in the results obtained with the MDPA and TSSP schemes. In both schemes, the simulated sea ice extent with concentration larger than 0.9 is smaller in winter as compared to results of the upwind scheme. The sea ice covered area in Bering Sea becomes smaller and the meridional gradient of the sea ice concentration in the margin of Okhotsk Sea becomes larger. In summer, the sea ice extent with sea ice concentration larger than 0.9 in the Greenland Sea and Chukchi Sea contracts northward. Although the performances of different advection schemes can be analyzed in simple cases, in real complicated applications, due to influences of other mechanisms and complicated interactions involved, more uncertainties may arise. References
Fig.7 As in Fig.1 except for IRS scheme
[1]
[2]
[3] [4]
Fig.8 As in Fig.2 except for IRS scheme
3. Conclusion In the study, we introduce three new advection schemes into the sea ice model SIS. With an introduction of a simple mixed layer ocean model into the SIS to calculate the variations of the sea surface temperature, four numerical simulations of the sea ice using different advection schemes are performed. Simulation results in the northern high latitudes are analyzed. It is shown that the SIS can reproduce the seasonal variations of the sea ice well with whichever advection scheme. In all simulations, the SIS produces some extra sea ice cover in the Greenland Sea and Barents Sea in winter due to the neglect of the influence of ocean currents. The sea ice model produces smaller ice thickness in the Arctic basin than the observation value, which may be due to an inappropriate parameterization of the heat flux from the lower ocean water. There are similar features among results obtained with the introduced three advection schemes. The simulated sea ice thicknesses are all smaller than that of the upwind scheme and the simulated sea ice
[5] [6] [7]
[8]
[9]
[10] [11] [12]
PARKINSON C. L., RIND D. and HEALY R. J. et al. The impact of sea ice concentration on climate model simulations with the GISS GCM[J]. J. Climate, 2001, 14(12): 2606-2623. ALEXANDER M. A., BHATT U. S. and WALSH J. E. et al. The atmospheric response to realistic Arctic sea ice anomalies in an AGCM during winter[J]. J. Climate, 2004, 17(5): 890-905. HOLLAND M., BIT C. and HUNKE E. et al. Influence of the sea ice thickness distribution on polar climate in CCSM3[J] J. Climate, 2006, 19(11): 2398-2414. SEDLÁCEK J., MYSAK L. A. Sensitivity of sea ice to wind-stress and radiative forcing since 1500: A model study of the little ice age and beyond[J]. Clim. Dyn., 2009, 32(6): 817-831. HOLLAND M. M., BITZ C. M. Polar amplification of climate change in coupled models[J]. Clim. Dyn., 2003, 21(3-4): 221-232 HILMER M., LEMKE P. On the decrease of arctic sea ice volume[J]. Geophys. Res. Let., 2000, 27(22): 37513754. WANG Gang, JI Shun-ying and LV He-xiang et al. Drucker-prager yield criteria in viscoelastic-plastic constitutive model for the study of sea ice dynamics[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(6): 714722. WANG H., SKAMAROCK W. C. and FEINGOLD G. Evaluation of scalar advection schemes in the advanced research WRF model using large-eddy simulations of aerosol-cloud interactions[J]. Mon. Wea. Rev., 2009, 137(8): 2547-2558. YANG Wei, LIU Shu-hong and WU Yu-lin. An unsplit lagrangian advection scheme for volume of fluid method[J]. Journal of Hydrodynamics, 2010, 22(1): 73-80. MARTIN T. Arctic sea ice dynamics: Drift and ridging in numerical models and observations[D]. Ph. D. Thesis, Potsdam, Germany: Alfred Wegener Institute, 2007. WINTON M. A reformulated three-layer sea ice model[J]. Journal of Atmospheric and Oceanic Technology, 2000, 17(4): 525-531. DELWORTH T. L. and Couuthors. GFDL’s CM2 global coupled climate models. Part I: Formulation and
378
[13]
[14]
[15] [16]
simulation characteristics[J]. J. Climate, 2006, 19(5): 643-674. SMOLARKIEWICZ P. K. Multidimensional positive definite advection transport algorithm: An overview[J]. International Journal for Numerical Methods in Fluids, 2006, 50(10): 1123-1144. LIU Xi-ying, ZHANG Xue-hong and YU Ru-cong et al. Fine-resolution simulation of surface current and sea ice in the Arctic Mediterranean Seas[J]. Chinese Journal of Oceanology and Limnology, 2007, 25(2): 132-138. DUKOWICZ J. K., BAUMGARDNER J. R. Incremental remapping as a transport/advection algorithm[J]. J. Comput. Phys., 2000, 160(1): 318-335. LIPSCOMB W. H., HUNKE E. C. Modeling sea ice transport using incremental remapping[J]. Mon.Wea. Rev., 2004, 132(6): 1341-1354.
[17] [18] [19]
[20]
STEELE M., MORLEY R. and ERMOLD R. PHC: A global ocean hydrography with a high quality Arctic Ocean[J]. J. Climate, 2001, 14(9): 2079-2087. UPPALA S. M., KÅLLBERG P. W. and SIMMONS A. J. et al. The ERA-40 re-analysis[J]. Quart. J. R. Meteorol. Soc., 2005, 131(612): 2961-3012. LIU X. Y., ZHANG X. H. and YU Y. Q. et al. Mean climatic characteristics in high northern latitudes in an ocean-sea ice-atmosphere coupled model[J]. Adv. Atmos. Sci., 2004, 21(2): 236-244. BOURKE R. H., GARRETT R. P. Sea-ice thickness distribution in the Arctic Ocean[J]. Cold Reg. Sci. Technol., 1987, 13(3): 259-280.