Implementation of a One-Dimensional Enthalpy Sea-Ice Model in a Simple Pycnocline Prediction Model for Sea-Ice Data Assimilation Studies

Xinrong WU,Shaoqing ZHANG,and Zhengyu LIU

1Key Laboratory of Marine Environmental Information Technology,State Oceanic Administration, National Marine Data and Information Service,Tianjin 300171

2Geophysical Fluid Dynamics Laboratory,National Oceanic and Atmospheric Administration, Princeton University,Princeton,NJ 08542,USA

3Center for Climate Research and Dept.Atmospheric and Oceanic Sciences, University of Wisconsin-Madison,Madison,WS 53706,USA

4Laboratory of Ocean-Atmosphere Studies,Peking University,Beijing 100871

Implementation of a One-Dimensional Enthalpy Sea-Ice Model in a Simple Pycnocline Prediction Model for Sea-Ice Data Assimilation Studies

Xinrong WU?1,Shaoqing ZHANG2,and Zhengyu LIU3,4

1Key Laboratory of Marine Environmental Information Technology,State Oceanic Administration, National Marine Data and Information Service,Tianjin 300171

2Geophysical Fluid Dynamics Laboratory,National Oceanic and Atmospheric Administration, Princeton University,Princeton,NJ 08542,USA

3Center for Climate Research and Dept.Atmospheric and Oceanic Sciences, University of Wisconsin-Madison,Madison,WS 53706,USA

4Laboratory of Ocean-Atmosphere Studies,Peking University,Beijing 100871

To further explore enthalpy-based sea-ice assimilation,a one-dimensional(1D)enthalpy sea-ice model is implemented into a simple pycnocline prediction model.The 1D enthalpy sea-ice model includes the physical processes such as brine expulsion,flushing,and salt diffusion.After being coupled with the atmosphere and ocean components,the enthalpy sea-ice model can be integrated stably and serves as an important modulator of model variability.Results from a twin experiment show that the sea-ice data assimilation in the enthalpy space can produce smaller root-mean-square errors of model variables than the traditional scheme that assimilates the observations of ice concentration,especially for slow-varying states.This study provides some insights into the improvement of sea-ice data assimilation in a coupled general circulation model.

sea ice,enthalpy,coupled model,data assimilation,ensemble Kalman filter

1.Introduction

Sea ice interacts with other components in the polar climate system(e.g.,Budyko,1969)and plays crucial roles in global climate(e.g.,Walsh,1983).Due to model error,the modeled sea-ice seasonal cycle and anomaly variability are different from observations.Because of the interactions of sea ice and the ocean-atmosphere system,the modeled seaice bias can contribute to coupled-model errors in both climate features and variability.One usually attempts to assimilate observations of sea ice into a coupled model to correct sea-ice states,so as to improve the quality of climate assessment and predictions.Since ice concentration is the most common observation of sea ice,most present sea-ice data assimilations are based upon this variable(e.g.,Lis?ter et al.,2003;Lindsay and Zhang,2006;Stark et al.,2008). However,there are two outstanding issues in traditional(i.e., present)sea-ice data assimilation methodologies.First,ice concentration is bounded to have values between 0 and 1,implying that its errors cannot follow a Gaussian distribution. This conflicts with the commonly held assumption in most data assimilation schemes that the errors in the background state and observations are Gaussian.In the ensemble Kalman filter(EnKF)context,when the sea-ice temperature is sufficiently low or high,all ensemble members of the ice concentration may converge to 1 or 0,causing abnormal error covariances between ice concentration and observations.Second,the assimilation update of ice states may be overridden by the strong impact of the atmospheric and oceanic forcings (e.g.,Toyoda et al.,2011).In order to overcome these difficulties of sea-ice data assimilation,Zhang et al.(2013)developedanenthalpy-basedsea-icedataassimilationschemewith a simplified coupled system that uses a nonlinear function as a high-order approximation of the sea-ice model.In the real world,however,sea ice involves many complex physical processes that have a substantial effect on the climate system. For example,the processes of brine expulsion and flushing modulate the salinity structure,which influences ocean circulation.Here,we introduce a more realistic one-dimensional (1D)enthalpy sea-ice model(Notz,2005;Notz and Worster, 2006)into a 1D atmosphere-ocean coupled model(Zhang,2011a,b).The sea-ice model employed here simulates the enthalpy and bulk salinity evolution of sea ice.The aim of the 1D coupled model developed in this paper is to enable the study of error-continuous sea-ice data assimilation.

?Institute of Atmospheric Physics/Chinese Academy of Sciences,and Science Press and Springer-Verlag Berlin Heidelberg 2016

The paper is structured as follows:Section 2 introduces the simple coupled model,and section 3 presents the numerical implementation and evaluation of the enthalpy sea-ice model.We examine the impact of sea ice on model variability in section 4,and present a case study of sea-ice data assimilation in section 5.A summary and further discussion conclude the paper in section 6.

2.The simple coupled model

The model used in this study is based on the simple coupled model developed by Zhang et al.(2013),in which a simple nonlinear sea-ice model is coupled with Lorenz's chaotic atmosphere(Lorenz,1963)and a slab upper ocean that interacts with the deep ocean(Gnanadesikan,1999).We replace the simple sea-ice model with a more realistic enthalpy seaice model(Notz and Worster,2006)in this study.The governing equations of the coupled model are

where x1,x2and x3are high-frequency variables,with the originalσ,κand b parameters representing the atmosphere, and w represents the anomaly of the upper ocean temperature.A dot above the variable denotes the time tendency. The imposed external forcing of the slab ocean is parameterized by Sm+Sscos(2πt/Spd),in which Smand Ssdefine the magnitudes of the annual mean and seasonal cycle of the forcings.The term(Om,Od)defines the oceanic time scale,-c7?v,t-1simulates the reduced radiative heat for the ocean due to the insulation of sea ice,andηrepresents the anomaly of deep-ocean pycnocline depth;its tendency equation is derived from the two-term balance model of the zonaltime mean pycnocline(Gnanadesikan,1999).The function Φ representstheforwardoperatoroftheenthalpysea-icemodel. The input arguments of Φ include the atmospheric(c8x2)and oceanic forcings(c9w),as well as the previous sea-ice state variables(?v,T,H and Sbu).Here,?v,T,H and Sburepresent the solid volume fraction,temperature,enthalpy per unit volume,and bulk salinity of sea ice,respectively.Among these four quantities,H and Sbuare prognostic variables,while?vand T are diagnostic variables.Note that in this simple system,the sea ice impacts the atmosphere(x2)only through the ocean-atmosphere(w-x2)interaction.The terms c8x2and c9w serve as the top and bottom anomaly forcings of ice temperature.Note that to simulate the seasonal variation of the sea ice,the seasonal solar forcings Sscos(2πt/Spd)and 0.1Sscos(2πt/Spd)are also imposed on the top and bottom boundary conditions of ice temperature.

A leap-frog time differencing scheme with a Robert-Asselin time filter(Robert,1969;Asselin,1972)is applied to the atmosphere and the upper-and deep-ocean components, while an explicit Euler forward scheme is used in the sea-ice model.The time step of the atmosphere and the ocean is denoted as ΔtAO,while the sea-ice time step is denoted as ΔtI.For simplicity,we use“time step”to represent ΔtAOin this study;and one time unit(denoted as TU)is equal to 100 time steps.Note that the unit of measurement of variable t in Eq.(1)is TU.We will determine the value of ΔtIin section 3. The coupled model is initialized with the parameter values as (σ,κ,b,c1,c2,c3,c4,c5,c6,c7,c8,c9,Om,Od,Sm,Ss,Spd,Γ,γ) =(9.95,28,8/3,0.1,1,0.01,0.01,1,0.01,0.1,0.025,0.01,10, 1,10,1,10,100,0.25)(Zhang et al.,2013),whereγis the time-filtering coefficient in the Robert-Asselin time filter.

3.Numerical implementation of the sea-ice model

With the enthalpy sea-ice model,the methodology of enthalpy-based sea-ice data assimilation can be realized explicitly.Inthissection,weimplementthe1Denthalpysea-ice model based on the prototype of Notz(2005).

Based on the governing equations(see Appendix A)and the desalination processes(see Appendix B),we first outline the configuration and integration of the 1D enthalpy sea-ice model,then examine the roles of the physical processes in the sea-ice evolution.Lastly,we briefly introduce the couplings between the sea ice and other model components.

3.1.Configuration and integration of the sea-ice model

For the first trial,the vertical depth is set to 2 m,which is evenly partitioned into 200 layers,and the thickness Δz of each layer is therefore 1 cm.Table 1 lists all the seaice model variables and their units,while Table 2 gives the values and meanings of the parameters in the seaice model.The maximum stable time step is given by theCourant-Friedrich-Levy criterion,which is formulated as 0.05(Δz)2ρlcl/max(ks,kl)=9.81 s.According to the parameter configuration in section 2,the period of the external forcing of the slab ocean(i.e.,the value of Spd)is 10 TUs.That is,10 TUs equal one model calendar year.If we roughly assume one year has 360 days,1 TU will be equivalent to 3.1104×106s.Thus,ΔtAOequals 31 104 s.To facilitate the couplingbetweenthesea-iceandothercomponentsandmake ΔtIless than 9.81 s,the sea-ice time step ΔtIis set to 8.64 s, which means that the sea-ice model should be integrated by 3600 steps once the atmosphere and ocean have been integrated by one step.Appendix C lists the sequential steps to integrate the sea-ice model from time t to time t+1,while Fig.1 shows the corresponding flow chart.

Table 1.List of sea-ice model variables.

Table 2.Default values and meanings of the parameters in the seaice model.

3.2.Examination of the roles of physical processes in the sea-ice evolution

3.2.1.Brine expulsion

Brine expulsion is the most important physical process in the freezing period.Here,we provide a schematic illustration of how brine expulsion affects the sea-ice state,especially the bulk salinity.Figure 2 shows an example of the brine expulsion process with one day's integration from the initial temperature(0?C)and bulk salinity(10 psu)profile without coupling the sea-ice model with other components,where the dashed and solid curves represent the initial and final states of T(Fig.2a),Sbu(Fig.2b),?v(Fig.2c)and H(Fig.2d). The surface forcing is initialized as-10?C and subsequently varied as a cosine function from-10?C to 10?C with a oneday period,simulating the diurnal cycle.After one day's integration,the surface salt is expelled downwards and the thermal variables(temperature and enthalpy)are adjusted accordingly.The ice thickness reaches about 6 cm.Note that the Sbuvalues less than 10 psu in Fig.2b are caused by the varying external forcing.At the beginning of the day,the surface forcing is lower than the freezing point,causing the occurrence of brine expulsion.During the middle part of the day,the surface forcing is higher than the freezing point,and thus the brine expulsion is not activated.At the end of the day,the surface forcing is once again lower than the freezing point,and thus brine expulsion happens again.

3.2.2.Flushing

Although flushing is not a simple 1D vertical process (Eicken et al.,2002;2004),we still include it into the enthalpy sea-ice model so as to simulate the evolution of bulk salinity in summer.It should be noted that since the time scale of the flushing in the 1D enthalpy sea-ice model is greatly different to those in other model components,the magnitude of the flushing speed of the melt water is artificially multiplied by 10-7to sustain the coupled system steadily.

Fig.1.Flow chart of the numerical implementation of the enthalpy sea-ice model for a forward step(time t to t+1).

Fig.2.One-day integration results(solid line)of(a)temperature,(b)bulk salinity,(c)solid volume fraction and (d)enthalpy per volume,demonstrating the role of the brine expulsion process.The dashed lines represent the initial sea-ice state.The surface forcing is initialized as-10?C and subsequently varied as a cosine function from-10?C to 10?C with a one-day period.

In order to evaluate the role of flushing in simulating the evolution of salinity in the sea ice,the variation of the surface Sbuwith a 3-yr integration of the sea ice is shown in Fig.3,where the solid and dashed curves represent the results without and with flushing,respectively.The surface forcing is initialized by-10?C and subsequently varied as a cosine function from-10?C to 10?C with a 1-yr period in this long run for demonstrating the effect of flushing.It is evident that without flushing,the surface Sbuis artificially increased when the melt pond stands over the sea ice.Once flushing is introduced,although it is idealized,the surface Sbucan be partially reduced until the surface is refrozen or no ice exists underneath.

3.2.3.Salt diffusion

In this study,salt diffusion is mainly used to smooth the salinity profile when two ice sheets hit each other during the refreezing time.Three-point smoothing with the weighted coefficientα[Eq.(B8)]is implemented.Figure 4 compares the Sbuprofiles after a 1-yr integration with the same initial fieldsandsurfaceforcingasthoseinFig.3fordifferentαvalues.Without salt diffusion,large oscillations caused by the interaction of the upper and lower pieces of ice exist.With large smoothing coefficients(e.g.,10-5,5×10-5and 10-4), the salinity profile is gradually smoothed.We choose an intermediate smoothing coefficient,5×10-5,to implement the salt diffusion in the 1D enthalpy sea-ice model in the following analyses of the coupled model.

3.3.The couplings between the sea-ice and other model components

The top boundary condition of the sea ice is a seasonal solar forcing[Sscos(2πt/Spd)]imposed on the atmospheric temperature(c8x2).The bottom boundary condition is set by a small seasonal solar forcing[0.1Sscos(2πt/Spd)]imposed on the upper ocean temperature(c9w).The solar forcing is added to reflect the seasonal variation of the sea ice,while c8x2and c9w are used to induce the variability of the sea ice. The feedback from the sea ice to the ocean is simulated by the term-c7?t-1inthewequation,whichrepresentsthereduced radiative heat of the ocean due to the albedo of sea ice andvariation of each model component.The initial values of (x1,x2,x3,w,η)are set to(0,1,0,0,0).The initial T and Sbuof sea ice are set to 0?C and 10 psu,while other quantities are computed accordingly.The enthalpy sea-ice model component is activated after a spin-up period(500 TUs)of the atmosphere and the upper and deep ocean.The total simulation period is 35000 TUs,i.e.,3500 years.

Fig.3.Variations of surface bulk salinity with(dashed line)and without(solid line)flushing for a 3-yr integration of the sea-ice model.The surface forcing is initialized by-10?C and subsequently varied as a cosine function from-10?C to 10?C with a 1-yr period in this long run for demonstrating the effects of flushing.

Fig.4.Bulk salinity profiles after a 2-yr integration with the same initial conditions and the surface forcing as that in Fig.3, for different smoothing coefficients(black,blue,green and red for 0,10-5,5×10-5and 10-4,respectively).

Power spectrum analysis is used to examine the impact of the sea ice on the coupled model variability.Time series of anomalies of x2,w andηbetween 100 and 3000 years,with andwithoutthesea-ice component,areusedto implementthe power spectrum analysis(Fig.5).It is evident that the introduction of sea ice can enhance the power spectrum for ocean variables.For example,the largest power spectrum of w(η)is augmented from 25(120)years to 100(280)years.Thus,in this particular coupled model setting,sea ice can greatly affect the climate scale variability of the coupled system.Figure 6 presents the time series of the surface T(Fig.6a),surface Sbu(Fig.6b),bottom T(Fig.6c)and bottom Sbu(Fig.6d) of the sea-ice model starting from the 50-yr spin-up of other components.The sea-ice temperature responds to the boundary forcings very quickly,while the bulk salinity has a much longer time scale than the temperature.Due to the inclusion of high-frequency atmospheric information in the upper forcing,the surface temperature of sea ice also varies quickly,in addition to the seasonal signal.Likewise,the low-frequency upper ocean leads the smoothly varying bottom temperature of sea ice.Because salinity exchange occurs after heat exchange,the time scale of salinity(about 1000 years here)in the sea ice is longer than that of the ice temperature.Figure 7 shows the final ice profiles of T(Fig.7a),Sbu(Fig.7b)and ?v(Fig.7c)when a melt pond lies over the sea ice,which means that flushing is happening.It is worth mentioning that the current simple model setting has limitations and can only serve as a preliminary proof-of-concept study on the enthalpy sea-ice model and the impact of enthalpy-based sea-ice data assimilation on climate estimation.More discussion on the limitations of the simple model setting and further follow-up studies to address these limitations will be given in section 6. the atmosphere being only indirectly affected by the sea ice through the air-sea interaction in the simple system.Note that?v,t-1here is the vertical average of?vin the sea-ice model at time t-1.

4.Model variability

Withthedefaultvaluesofthemodelparametersdescribed in sections 2 and 3,the model is first integrated to study the

5.Enthalpy-based sea-ice data assimilation

In this section,with the newly developed 1D enthalpy sea-ice model,we study the impact of enthalpy-based seaice data assimilation on coupled state estimation.To roughly simulate the situation in the real world(climate models are always biased compared with observations),we use a biased twin experiment framework.

5.1.Biased twin experiment setup

Fig.5.Power spectra of the anomalies of(a,b)x2,(c,d)w,and(e,f)η,without(left panels)and with (right panels)the sea-ice component.Dashed curves represent the 95%upper confidence limits.Note that due to the short time scale of x2,the unit of measurement on the x-axis in(a,b)is time-step and one year is equal to 1000 time-steps.

The standard model setting described in sections 2 and 3.1 defines the“truth”model that is used to produce“observations”.Starting from the states at the end of the 3500-yr integration,the truth model is further integrated for another 100 years to generate the true solution of the assimilation problem.A Gaussian noise is superimposed on the“truth”to simulate the observational error.We assume that only x1,2,3(representing x1,x2and x3),w and the surface?vare observed,and the standard deviations(relative to truth) of these observations are 2.0(for x1,2,3),0.5(for w)and 0.05 (for?v).Note that if the observation of?vis greater than 1 or less than 0,it will be set to 1 or 0.The observational interval of the atmosphere is 5ΔtAO,while the upper ocean and seaice observational interval is 10ΔtAO.The biased assimilation model is defined by setting a differentκparameter value(29) from the standard model setting.Then,following the“truth”model spin-up procedure,the assimilation model is also integrated by 3500 years,starting from the same initial conditions.The ensemble initial conditions from which the ensemble filtering data assimilation starts are produced as follows: A Gaussian random noise with the same standard deviation of x1,2,3(w)observational error is added to the atmospheric (oceanic)states to form the ensemble initial condition of the atmosphere(ocean).The ice states are not perturbed.

Fig.6.Time series of(a)surface temperature,(b)surface bulk salinity,(c)bottom temperature and(d) bottom bulk salinity of the sea-ice model starting from the states at the end of the 500-TU spin-up without the sea-ice component.Note that due to the different time scales of ice temperature and bulk salinity, panels(b,d)use a longer time period than panels(a,c).

To understand the impact of enthalpy-based sea-ice data assimilation on coupled state estimation,we perform the following three data assimilation experiments:(1)a traditional coupled data assimilation experiment without sea-ice data assimilation(denoted as EXP1);(2)additional to(1),an experiment that assimilates the observations of surface?vand surface T to respectively adjust the modeled surface?vand surface T(denoted as EXP2);and(3)additional to(1),an experiment that assimilates the observations of surface H into the model(denoted as EXP3).For simplicity,we adopt the simplest assimilation scheme for the sea-ice assimilations in this study.In EXP1,observations of x1,2,3and w are used to adjust the modeled x1,2,3and w in a multivariate way.In EXP2 and EXP3,observations of x1,2,3and w are first used to adjust the modeled x1,2,3and w,as in EXP1.Then,the surface sea-ice observations are used to adjust surface ice states in a univariate way.Note that if a negative value of?v,or a value larger than 1,arises during the assimilation of?vin EXP2,it will be set to 0 or 1.All assimilation experiments adopt the same initial condition and observing system,and are conducted for 100 years.

Fig.7.Profiles of(a)temperature,(b)bulk salinity and(c)solid volume fraction at the end of the 3500-yr integration.

For the assimilation of surface H,when observations of surface?vand x2are available,the observation of surface T is first derived from x2according to the coupling from the atmosphere to the sea ice.Then,the observation of surface ?mis computed by Eq.(A16).Note that observed?mis also bounded between 0 and 1.Afterwards,the observation of surface H is computed using Eq.(A5).Finally,the observation of surface H is used to adjust the modeled surface H. To reasonably set the standard deviation of the observational error of H,we first calculate the natural variabilities(i.e.,the climatological standard deviation)of the surface?vand H as 0.27 and 8.6×107J m-3,respectively,using the last 1000-TU model integration of the truth model.Because the standard deviation of the surface?vobservational error is set to 0.05,the ratio of the observational error to the natural variability for the surface?vis 0.19.To achieve a similar ratio for the surface H,we set the standard deviation of the surface H observational error as 1.6×107J m-3.

According to the above description,EXP3 actually uses two types of observations:the observation of surface?v,and the derived observation of surface T.To conduct a fair comparison between EXP2 and EXP3,besides the assimilation of the surface?v,we also assimilate the derived observation of surface T in EXP2.

The data assimilation method adopted here is the ensemble adjustment Kalman filter(EAKF,Anderson,2001,2003), which is a derivative of the EnKF(Evensen,1994).Following previous studies(Zhang and Anderson,2003;Zhang, 2011a,b),theensemblesizeissetas20throughoutthisstudy. A two-time level adjustment(Zhang et al.,2004)is implemented in the state estimation.

5.2.Assimilation results

We first examine the probability density functions(PDFs) of the surface H and surface?vand display them in Fig.8. The H PDF in the 1D enthalpy sea-ice model shows a bimodal distribution,meaning the sea-ice model possess strong nonlinearities.The?vPDF is discontinuous at the points of ?v=0 and?v=1,while the H PDF is continuous(the end of the right-hand side is determined by the range of x2values when the temperature is above the freezing point).The discontinuous?vPDF can cause significant difficulties in correctly evaluating the associated error covariance in the data assimilation.Next,we compare the assimilation results of the three experiments.

Figure 9 presents the time series of the“truth”(red)and the analysis ensemble means of x2(Fig.9a),w(Fig.9b),η (Fig.9c),the surface H(Fig.9d),the surface Sbu(Fig.9e) and the surface?v(Fig.9f)in EXP1(black),EXP2(green) and EXP3(blue).On the one hand,compared to the results of EXP1,which does not use any sea-ice observation,EXP2 converselyworsensallmodelvariables.Due tothestrongimpact of the atmospheric forcing,some(or all)of the ensemble members of the surface?vmay be 0 or 1,causing the poorly evaluated error covariance between the modeled surface?vand the observed surface?v,which further contaminates the assimilation quality of EXP2.On the other hand,due to the continuous property of the H-associated error covariance,the enthalpy-based sea-ice data assimilation in EXP3 can significantly improve the quality of estimated model states.This is also evidenced by the measure of root-mean-squares error (RMSE)presented in Fig.10.

From Fig.10,we find that while the enthalpy-based seaice assimilation greatly improves the accuracies of the estimated surface H and surface?v,it also ameliorates the deepocean state.This could be attributable to the strong“memory”of the low-frequency nature of the deep ocean;otherwise,the deep-ocean variableηwould drift away due to the accumulation of errors.For EXP2,the results are very unstable,especially for the oceanic variables.It is worth men-tioning that,due to the strong nonlinearity of the 1D enthalpy sea-ice model,the performance of the sea-ice assimilation is alittledifferentfromtheidealizedcaseinZhangetal.(2013). First,the assimilation spin-up time is much longer.And second,it is not guaranteed that the state-estimation error in the enthalpy-based sea-ice assimilation is always the smallest.

Fig.8.The probability density functions(PDFs)of(a)surface enthalpy(108J m-3)and(b)surface concentration.The statistics are based on the last 1000-TU model integration of the“truth”model.

6.Summary and discussion

A 1D enthalpy sea-ice model is implemented into a simple climate model for further conducting enthalpy-based seaice data assimilation studies.The 1D enthalpy sea-ice model includes the physical processes of brine expulsion,flushing, and salt diffusion.By coupling with the atmosphere and ocean components,the enthalpy sea-ice model can be integrated stably and serves as an important modulator of model variability.Results from a biased twin experiment show that, due to the continuous PDF of enthalpy,the sea-ice data assimilation in enthalpy space can improve coupled state estimates.Theresultsofthisstudy providesomeinsightsintothe improvementofsea-icedataassimilationinacoupledgeneral circulation model(CGCM).

The current work with a simple coupled model setting can only serve as a preliminary study on the enthalpy seaice model and the impact of enthalpy-based sea-ice data assimilation on climate estimation.Many limitations exist,and further studies need to address these limitations and increase our understanding.First,although Lorenz's chaotic atmospheric model is simplified from the finite-amplitude atmospheric convection equations(e.g.,Saltzman,1962)and thus has some representation of the features of atmospheric flows, the extreme simplification sets a limited frequency band for the model(Lorenz,1984).Second,the simplification of the ocean with a slab upper ocean interacting with a singlevariable deep ocean is also a limitation in terms of its representation of the real climate system,which has a rich spectrum.Third,due to the highly chaotic nature of the Lorenz atmosphere,the variability of the sea-ice model is dominated by that of the chaotic atmosphere,which is different from the realworld.Fourth,the overlysmall flushing speedofthe melt water may lead to the overly long time scale of bulk salinity in sea ice.Fifth,the coupled model does not strictly obey the basic laws of physics,especially with respect to energy and mass conservation,which is far away from a CGCM. And lastly,the sea-ice concentration in this study is different from that traditionally being assimilated.The traditional seaice concentration denotes the areal fraction covered by ice in a certain area,such as a 2D horizontal model grid,which can be easily observed by satellite.Whereas,in this paper, sea-ice concentration denotes the solid ice mass/volume fraction in the vertical column,such as an ice core.This sea-ice concentration is generally very difficult to observe,and also has a distinct solar radiation effect.Thus,before the application of the enthalpy-based sea-ice data assimilation,a conversion from the traditional observed sea-ice concentration to the concentration in this study should be established.For example,a relationship between the surface ice concentration in this study and the satellite-observed ice concentration and thickness could be set up.

There are also many challenges in implementing an enthalpy sea-ice model into a CGCM.Currently,the gravity drainage process is not yet implemented.Given the importance of the gravity drainage process in the change of sea ice states(e.g.,Notz,2005),the implementation of this process is the first priority in follow-up studies,and the impact of gravity drainage on sea-ice behavior in a simple model framework must be examined.Then,the 1D model must be expanded to 3D and linked with a coupled ocean-atmosphere system by taking sea-ice dynamics(e.g.,Flato and Hibler III, 1990;Pritchard et al.,2010)into account.In addition,in the current simple model development and study,we use 200 ver-tical levels for the 2-m sea-ice layer.In a CGCM case,such a sea-ice model layout demands an impractical computational cost.How to design the ice model layout so that the enthalpy sea-ice model does not tremendously increase the computational cost is another big challenge for CGCM application. Furthermore,the enthalpy is non-observable and cannot be perfectly defined in observational space,as in a model.Under this circumstance,another alternative continuous quantity (like sea-ice temperature)can be considered to implement the sea-ice data assimilation.

Fig.9.Time series of the analysis ensemble means of the(a)atmosphere variable x2,(b)upper ocean w, (c)deep oceanη,(d)surface enthalpy of sea ice,(e)surface bulk salinity of sea ice,and(f)surface solid volume fraction of sea ice in EXP1(black),EXP2(green)and EXP3(blue).The red curve represents the“truth”.Note that due to the different time scales,each panel uses different horizontal tick marks.

Fig.10.Time series of the root-mean-square errors(RMSEs)of the analysis ensemble mean of the(a) upper ocean w,(b)deep oceanη,(c)surface concentration,and(d)surface enthalpy in EXP1(black line),EXP2(green line)and EXP3(blue line).Note that due to the different time scales,each panel uses different horizontal tick marks.

Acknowledgements.Special thanks go to Dr.Dirk NOTZ at Max Planck Institute for his generous help during the model implementation.This research is co-sponsored by grants from the National Natural Science Foundation(Grant Nos.41206178, 41306006,41376015,41376013 and 41176003),the National Basic Research Program(Grant No.2013CB430304),the National High-Tech R&D Program(Grant No.2013AA09A505),and the Global Change and Air-Sea Interaction Program(Grant No.GASI-01-01-12)of China.

APPENDIX A

Description of the 1D Enthalpy Sea-Ice Model Equations

In his doctorate thesis,Notz(2005)gives a detailed description of a 1D enthalpy sea-ice model that is vertically partitioned into N layers.The prognostic variables include enthalpy and bulk salinity,while other model variables,like temperature,brine salinity and solid fraction,are diagnostic variables.We review the definition of enthalpy and the forwarding equations in this appendix.

A1.Definition of enthalpy

Generally,there are two kinds of enthalpy;namely,the enthalpy per unit mass(Hm)and the volume of a mushy layer (Hv).We first give the definition of Hm,which is usually formulated as

where the indices s and l denote solid and liquid phases,respectively;?mis the solid mass fraction;and Hsand Hlrepresent the enthalpies per unit mass of pure solid and liquid solution,respectively.Note that?m,Hsand Hlare functions of the temperature(in?C)of the mushy layer(hereafter denoted by T).

If a sea-ice sample has no solid fraction,i.e.,when thetemperature T is larger than the freezing point Tf,Hmequals Hl,which is calculated by T×cl,where cldenotes the specific heat capacity of the liquid.

When T is smaller than Tf,after introducing the latent heat L(T)=Hl(T)-Hs(T)of fusion,Hmcan be written as

If the enthalpy of the pure liquid solvent at the freezing point Tfis set to zero,Hsat temperature T is given by(Notz,2005)

where L0is the latent heat of fusion at the freezing point.Setting the heat capacity of the solid ice(cs)to be the value(cs,0) at Tfand approximating L in Eq.(A2)by L0,Eq.(A2)can be further simplified as

For the enthalpy per unit volume of a mushy layer(Hv), it is computed by Hmˉρ,whereˉρ=(1-?v)ρl+?vρsdenotes the volume-fraction-weighted average of the densities of the solid ice(ρs)and the liquid brine(ρl),and?vrepresents the solid volume fraction.Note that the computational formulation ofˉρa(bǔ)nd the relationship between?vand?mwill be deduced later.

A2.Governing equations

A2.1.The enthalpy forward equation

The enthalpy per unit volume Hv(hereafter denoted as H for brevity)is determined by the heat balance equation(see Alexiades and Solomon,1993):

Here,Q is the heat flux at the internal grid cell or the boundary layer,which is calculated as

where k is the heat conductivity(also represents heat diffusion coefficient)of sea ice,and is calculated as the volumefraction-weighted average of heat conductivities of the solid ice(ks)and the liquid brine(kl):

The boundary heat fluxes from the atmosphere and the ocean act as the upper and the lower boundary conditions of Q,respectively.

A2.2.The salinity equations

There are two kinds of sea-ice salinity.One is the brine salinity Sbr,which denotes the salinity of the liquid brine,and is formulated as

where mldenotes the mass of the liquid brine and msaltrepresents the salt of a sea-ice sample.In sea-ice studies,Sbris usuallyobtainedfromobservationsofT withalinearliquidus relationship.An often used linear expression(e.g.,Notz, 2005)of Sbris as follows:

The other kind of sea-ice salinity is the bulk salinity Sbu, which represents the salinity of the water that is formed when a sea-ice sample is completely molten.Sbuis a common measure of the salinity in sea ice.The formula of Sbuis

where msdenotes the mass of the solid pure ice.As noted by Notz(2005),Sbuis usually much lower than that of the water from which it has formed and most of the salt is lost during the formation and growth of the ice.Thus,the evolution of Sbustrongly correlates to the desalination processes (Untersteiner,1968),such as brine expulsion,flushing,salt diffusion and gravity drainage.In this study,the latter process is not included in the enthalpy ice model.

To forecast Sbu,we first define the absolute salt content for each layer as Sabs=SbuΔzˉρ.Since the desalination processes only affect the salinity in the brine,thus the update of Sbuis realized after the update of Sabs,which is adjusted by Sbr.See Appendix B for the detail of the update of Sabs.

A2.3.The ice concentration equations

There are two key quantities of ice concentration:the solid mass fraction?mand the solid volume fraction?v.The former is defined as the solid mass per mass of mush,i.e.,

AccordingtoEqs.(A9)and(A11),Eq.(A12)canberewritten as

For?v,the definition leads to the following formula:

where Vsand V represent the volumes of the solid and the sea ice,respectively.Accordingly,ρs=ms/(?vV)andρl= ml/[(1-?v)V].Thus,the mean density can be formulated as

Finally,?mcan be written as

So,the expression of?vis

A2.4.The diagnosis of temperature

T is diagnosed by Hmand Sbuas follows:If Hmis larger than Tf×cl,T is directly diagnosed by Hm/cl.Otherwise, substituting Eq.(A10)and Eq.(A13)into Eq.(A4),with the knownbulksalinitySbuandHm,wehaveaquadraticequation of T as follows:

Using a Newton-Raphson method,T can be solved numerically.

APPENDIX B

Physical Processes in the 1D Enthalpy Sea-Ice Model

B1.Brine expulsion

Brine expulsion is the most important desalination process(Cox and Weeks,1975)and is driven by the pressure gradient between brine cavities caused by freezing and melting processes.The velocity field UUU induced by the brine expulsion can be derived from the following mass-conservation law:

With the expression ofρˉand the approximation ofρland ρsbeing constants,the divergence of UUU can be written as (Worster,1992)

In the discrete formation,the net displacement of the brine for the ith mushy layer during one time-step can be written as

where Δzidenotes the thickness of the ith layer.

B2.Flushing

Flushing happens when melt water stays over ice,which drives the brine downward.This mechanism occurs only in summer,which contributes to the salinity profile of multiyear sea ice(Untersteiner,1968).The velocity of the flushing can be calculated from Darcy's law(e.g.,Eicken et al.,2004):

whereΔzmprepresentsthedownwarddisplacementofthesurface of the melt pond,μis the dynamic viscosity of the brine, g is gravity,ΔtIis the time-step of the sea-ice model,h is the ice thickness,Πminis the minimum specific permeability of the mushy layer,which is formulated as

(Freitag,1999),and zmpis the hydraulic head,which is computed as(Notz,2005)

where“front”represents the interface layer between the ice and the sea water.In discrete formation,the displacement of the brine caused by the flushing for the ith layer is simply given by

which together with Eq.(B3)constitutes the net displacement of the ith layer.

B3.Salt diffusion

Due to the existence of gravity drainage,the liquid phases (i.e.,the brines)within a mushy layer interact with each other through salt diffusion.Thus,similar to brine expulsion and flushing,salt diffusion affects bulk salinity through the absolute salt content.Although the related theory(e.g.,Mellor and Kantha,1989;Notz,2005)of salt diffusion exists,to ease the complexity of the ice model and not lose the generality, we simply use a smoothing scheme to represent the salt diffusion process,which is especially designed for the situation when two pieces of ice hit each other at a certain layer.The formulation is

where i is the index for the mushy layer,which varies from 2 to N-1,andαrepresents the smoothing coefficient.

APPENDIX C

Sequential Steps of the Implementation of the 1D Enthalpy Sea-Ice Model

Step 1:After defining the parameters,the T(temperature)and Sbu(bulk salinity)fields are initialized.Then,the freezing-point(Tf)is computed using Tf=-mSbu,where m represents the solid mass.The liquidus relationship[Eq. (A10)]is used to calculate the brine salinity(Sbr),which is combined with Sbuto initialize the solid mass(?m)and volume fraction(?v)with Eqs.(A13)and(A17).Afterwards, the absolute salt content(Sabs)and the enthalpy per unit mass(Hm)are initialized.Lastly,the heat conductivity(k)is initialized with Eq.(A8).

Step 2:Compute the heat flux(Q)with Eq.(A7),which includes the upper and lower temperature forcings.

Step 3:Update the enthalpy per unit volume(H)by integrating Eq.(A6),after which the absolute enthalpy content (Habs)of a certain grid cell is computed using Habs,i=ΔziHi.

Step 4:For a given vertical brine displacementδzifor the ith layer,the absolute enthalpy content of this grid cell (Habs,i)is adjusted by

where Hl,i=clρlTiis the enthalpy of the moving brine in the ith layer,and j is either+1 or-1,depending on whether the flow through a cell is upward or downward.Then,H is updated again with Hi=Habs,i/Δzi.

Step 5:With T,?vand Sbuat the last time step,Sbris updated using the liquidus relationship when?vis greater than 0;otherwise,it equals to Sbu.

Step 6:Similar to Step 4,the absolute salt content for the ith layer(Sabs,i=Sbu,iΔziˉρi)changes through import or export of salt with moving brine according to(Notz,2005)

Step 7:If the top grid cell becomes fully liquid(?v,1=0) but a certain grid underneath has sea ice[(max(?v)＞0)], flushing occurs.The displacement of each grid cell is computed with Eq.(B7).Then,the thickness of the first layer is reduced correspondingly.If it is less than 0,Habsand Sabsfields are shifted upwards by the size of the top grid cell,and H is also updated afterwards.

Step 8:Diagnose Sbuwith Sbu,i=Sabs,i/(Δziˉρi),which is used to compute the new Tf.If H≥Tfclρl,T is directly diagnosed as H/(clρl).Otherwise,compute Hmas H/ˉρ,then retrieve T from Eq.(A18)with the known Sbuand Hm.

Step 9:Perform salt diffusion according to Eq.(B8).

Step 10:With the new Sbuand Sbr,update?mand?vby Eq.(A13)and Eq.(A17).

Step 11:Update k for each layer with the updated?vby Eq.(A8).

Step 12:Compute the net displacementδzicaused by the brine expulsion and the flushing.

Note that Step 1 is the initialization,while Step-2 to Step-12 are the sequential procedures from time t to time t+1.

REFERENCES

Anderson,J.L.,2001:An ensemble adjustment Kalman filter for data assimilation.Mon.Wea.Rev.,129,2884-2903.

Anderson,J.L.,2003:A local least squares framework for ensemble filtering.Mon.Wea.Rev.,131,634-642.

Alexiades,V.,and A.D.Solomon,1993:Mathematical Modelling of Melting and Freezing Processes.Taylor&Francis,323 pp.

Asselin,R.,1972:Frequency filter for time integrations.Mon. Wea.Rev.,100,487-490.

Budyko,M.I.,1969:The effect of solar radiation variations on the climate of the Earth.Tellus,21,611-619.

Cox,G.F.N.,and W.F.Weeks,1975:Brine drainage and initial salt entrapment in sodium chloride ice.CRREL Res.Rept. 345,U.S.Army Cold Reg.Res.and Eng.Lab.,Hanover,NH, US.

Eicken,H.,H.R.Krouse,D.Kadko,and D.K.Perovich.,2002: Tracer studies of pathways and rates of meltwater transport through Arctic summer sea ice.J.Geophys.Res.,107(C10), SHE 22-1-SHE 22-20.doi:10.1029/2000JC000583.

Eicken,H.,T.C.Grenfell,D.K.Perovich,J.A.Richter-Menge, andK.Frey.,2004:HydrauliccontrolsofsummerArcticpack ice albedo.J.Geophys.Res.,109(C8),C08007,doi:10.1029/ 2003JC001989.

Evensen,G.,1994:Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.J.Geophys.Res.,99,10 143-10 162.

Flato,G.M.,and W.D.Hibler III,1990:On a simple sea-ice dynamics model for climate studies.Annals of Glaciology,14, 72-77.

Freitag,J.,1999:Untersuchungen zur Hydrologie des arktischen Meereises:Konsequenzen f¨ur den kleinskaligen Stofftransport=The hydraulic properties of Arctic sea-ice:Implications for the small scale particle transport,Berichte zur Polarforschung(Reports on Polar Research),Bremerhaven,Alfred Wegener Institute for Polar and Marine Research,No.325, 150 pp.

Gnanadesikan,A.,1999:A simple predictive model for the structure of the oceanic pycnocline.Science,283,2077-2079.

Lindsay,R.W.,and J.Zhang,2006:Assimilation of ice concentration in an ice-ocean model.J.Atmos.Oceanic Technol.,23, 742-749.

Lis?ter,K.A.,J.Rosanova,and G.Evensen,2003:Assimilation of ice concentration in a coupled ice-ocean model,using the ensemble Kalman filter.Ocean Dyn.,53,368-388.

Lorenz,E.N.,1963:Deterministic nonperiodic flow.J.Atmos. Sci.,20,130-141.

Lorenz,E.N.,1984:A very narrow spectral band.Journal of Statistical Physics,36,1-14.

Mellor,G.L.,and L.H.Kantha,1989:An ice-ocean coupled model.J.Geophys.Res.,94(C4),10 937-10 954.

Notz,D.,2005:Thermodynamics and fluid-dynamical processes in sea ice.PhD dissertation,University of Cambridge,238 pp.

Notz,D.,and M.G.Worster,2006:A one-dimensional enthalpy model of sea ice.Annals of Glaciology,44,123-128.

Pritchard,R.S.,G.Li,and G.F.N.Cox,2010:A simple sea ice drift forecast model.20th IAHR International Symposium on Ice,Lahti,Finland,12 pp.

Robert,A.,1969:The integration of a spectral model of the atmosphere by the implicit method.Proc.WMO/IUGG Symposium on NWP.Japan Meteorological Society,Tokyo,Japan,19-24.

Saltzman,B.,1962:Finite amplitude free convection as an initial value problem-I.J.Atmos.Sci.,19,329-341.

Stark,J.D.,J.Ridley,M.Martin,and A.Hines,2008:Sea ice concentration and motion assimilation in a sea ice-ocean model. J.Geophys.Res.,113,C05S91,doi:10.1029/2007JC004224.

Toyoda,T.,and Coauthors,2011:Impact of the assimilation of sea ice concentration data on an Atmosphere-Ocean-sea ice coupled simulation of the Arctic Ocean climate.Scientific Online Letters on the Atmosphere(SOLA),7,37-40,doi: 10.2151/sola.2011-010,037-040.

Untersteiner,N.1968:Natural desalination and equilibrium salinity profile of perennial sea ice.J.Geophys.Res.,73(4),1251-1257.

Walsh,J.E.,1983:The role of sea-ice in climate variability:Theories and evidence.Atmosphere-Ocean,3,229-242.

Worster,M.G.,1992:The dynamics of mushy layers.Interactive dynamics of convection and solidification:proceedings of the NATO Advanced Study Institute held at Chamonix,S.H. Davis et al.,Eds.,France.London,Kluwer Academic,113-138.

Zhang,S.,and J.L.Anderson,2003:Impact of spatially and temporally varying estimates of error covariance on assimilation in a simple atmospheric model.Tellus A,55,126-147.

Zhang,S.,J.L.Anderson,A.Rosati,M.J.Harrison,S.P.Khare, and A.Wittenberg,2004:Multiple time level adjustment for data assimilation.Tellus A,56,2-15.

Zhang,S.,2011a:Impact of Observation-optimized model parameters on decadal predictions:Simulation with a simple pycnocline prediction model.Geophys.Res.Lett.,38,L02702,doi: 10.1029/2010GL046133.

Zhang,S.,2011b:A study of impacts of coupled model initial shocks and State-Parameter optimization on climate predictions using a simple pycnocline prediction model.J.Climate, 24,6210-6226,doi:10.1175/JCLI-D-10-05003.

Zhang,S.,M.Winton,A.Rosati,T.Delworth,and B.Huang, 2013:Impact of enthalpy-based ensemble filtering sea-ice data assimilation on decadal predictions:Simulation with a simple pycnocline prediction model.J.Climate,26,2368-2378,doi:10.1175/JCLI-D-11-00714.1.

Wu,X.R.,S.Q.Zhang,and Z.Y.Liu,2016:Implementation of a one-dimensional enthalpy sea-ice model in a simple pycnocline prediction model for sea-ice data assimilation studies.Adv.Atmos.Sci.,33(2),193-207,

10.1007/s00376-015-5099-2.

17 April 2015;revised 13 July 2015;accepted 31 July 2015)

?Xinrong WU

Email:xinrong wu@yahoo.com