A coupled ice–mixed layer–ocean model is constructed for the Arctic Ocean, the Barents Sea, and the Greenland–Iceland–Norwegian Sea. The model is used to address Arctic numerical modeling with and without climate restoring. The model without climate restoring reproduces basic observed features of the Arctic ice–ocean circulation. The simulated oceanic processes adjust to the surface and lateral fluxes and transport heat and mass in a way that achieves a rough salt and heat balance in the Arctic in the integration period of seven decades. The main deficiency of the model is its prediction of unrealistically high salinity in the central Arctic, which tends to weaken the ocean currents.
The introduction of corrective salinity and temperature restoring terms has a significant impact on prediction of the ice–ocean circulation in the Arctic. The impact results from a chain reaction. First, the restoring terms change the salinity and thermal states in the oceanic mixed layer and below. The altered density structure, in turn, influences the ocean circulation by altering the ocean’s dynamic and thermodynamic processes. The ocean circulation then affects ice motion and ice thickness by altering the dynamics and thermodynamics of the ice. Restoring only ocean salinity induces a heat surplus or deficit, which causes the oceanic thermal state to drift away from the climatology. This is also the case with restoring both salinity and temperature in the upper ocean. Only restoring salinity and temperature in the deeper ocean is likely to avoid climate drift in the Arctic.
This paper points out the problems of using models without climate restoring and the consequences of using models with climate restoring for Arctic ice–ocean modeling. It stresses the need to use a realistic representation of surface fluxes in order to improve prognostic simulations. It also stresses the importance of oceanic processes such as convective overturning and vertical advection both in the Arctic and in the adjacent oceans. If these processes are not reasonably represented, the models can predict overwarm intermediate layers in the Arctic and an excessive heat influx at Fram Strait, causing the Arctic ice pack to possibly shrink gradually, a scenario global climate models should avoid.
The Arctic Ocean is thought to play a significant role in the changing earth climate system. A prominent feature of the Arctic is the presence of an ice cover over large areas. The ice cover substantially alters the transfer of heat, salt, and momentum between the atmosphere and the ocean and hence has the potential to alter atmospheric and oceanic circulation. The oceanic circulation, in turn, affects the ice dynamics and thermodynamics and the ice’s seasonal growth and decay. Because there are many uncertainties about physical processes (such as the surface heat, salt/freshwater, and momentum fluxes) at the interface between the ice and the atmosphere or the ocean and the atmosphere, numerical studies of the Arctic ice–ocean system often introduce climate-restoring terms (or damping terms, as they are often called) into the ocean models to correct errors in surface fluxes and other errors that the models may induce.
In fact, this approach is not limited to the Arctic Ocean modeling community. It is rather common with general circulation ocean models, which often use the restoring procedure one way or another. A number of global-scale ocean model studies (e.g., that of Sarmiento and Bryan 1982), use “robust” diagnostic ocean models that constrain ocean salinity and temperature to Levitus (1982) climatological data with a several-year restoring time. This loose restoring time allows the ocean model to “assimilate” observed data as well as simulate seasonal and short-timescale variability, while preventing the model from drifting away from ocean climatology on longer timescales. Meanwhile some global climate models constrain the ocean salinity and temperature by using flux correction at the interface of the atmosphere and the ocean to avoid climate “drift” (e.g., Manabe and Stouffer 1988; Sausen et al. 1988; Power 1995).
This restoring procedure was employed to varying degrees by Hibler and Bryan (1987), whose coupled ice–ocean model of the Arctic constrains ocean salinity and temperature to Levitus climatological data below a mixed layer of fixed depth. This restoring scheme was adopted in some later studies (e.g., Ries and Hibler 1991; Riedlinger and Preller 1991; Piacsek et al. 1991; Hibler and Zhang 1993; Hibler and Zhang 1994). Other studies, such as the work by Holland et al. (1996) and W. Maslowski (1997, personal communication), constrain only the surface values of salinity and temperature. The study by Piacsek et al. (1991) constrains both the surface salinity and temperature to climatology in addition to using climate restoring in the deeper ocean. One more study worth mentioning is that by Aukrust and Oberhuber (1995), which constrains only the surface salinity to climatology. Some details about the restoring schemes used in these studies are shown in Table 1.
Most of these restoring schemes rather effectively constrain the simulated ocean to the climatological mean state, and the inclusion of a diagnostic oceanic component generally improves the calculation of ice thermodynamics. However, the simulated ocean has a heavy bias toward the historic state, which poses questions about the state of the Arctic climate under conditions different from those in history and at present. Furthermore, the Arctic climatology is known to have a bias in favor of summer data, which may result in an excessively fresh upper ocean.
Consequently, a number of previous studies have attempted to use so-called prognostic ice–ocean models for the Arctic in which the climate-restoring terms are removed. Semtner (1987), for example, carried out a 20-yr integration of the Arctic ice–ocean system using an ice model with a simplified dynamics coupled to an ocean model with a fixed mixed layer. This same model was further used (Fleming and Semtner 1991) to show the importance of including a prognostic ocean in seasonal ice predictions. Meanwhile, a study by Ranelli and Hibler (1991) using a prognostic ice–ocean model of the Arctic with a 30-yr integration showed the importance of surface salt flux in the Arctic Ocean circulation and amplified the importance of better vertical resolution in shallow regions and the need for treating penetrative convection. More recently, a prognostic ice–ocean model was used by Hakkinen (1993), also for a 20-yr integration, to investigate the Arctic salinity anomaly. Hakkinen’s study showed a salinity drift of 0.1 ppt/yr in the Canadian Basin, which may also indicate the importance of using a realistic surface salt flux and a longer timescale for a prognostic model to reach a steady state. Note that the models used in most of these previous prognostic studies (with the exception of Hakkinen 1993) have relatively crude resolution (110–80 km). Resolution within this range does not well resolve the complex topography, coastlines, and narrow currents and passages of the ocean; this may not pose a serious problem in a short-term diagnostic study but can be problematic for a prognostic study, which often needs a sustained simulation over a longer timescale. The studies of Ranelli and Hibler (1991), Ranelli (1991), and Hakkinen (1993), for example, show that if a model does not signal a strong presence of Atlantic water in the Arctic, the density structure of the Arctic may eventually be changed so as to derive an unrealistically drifted ice–ocean system.
In this study, we constructed a higher-resolution, coupled, prognostic ice–mixed layer–ocean model for the ice–ocean system of the Arctic, Barents, and Greenland–Iceland–Norwegian (GIN) Seas and examined the characteristics of the model based on a 70-yr integration. It is a prognostic model in the sense that the ocean’s temperature and salinity are not constrained to climatological values and are forced only at the open boundaries. Meanwhile, we also carried out a number of sensitivity studies using similar models with different restoring schemes. Most of the restoring schemes are similar to some of those listed in Table 1. These diagnostic models, together with the standard prognostic model, were run repeatedly for ten times over the 7-yr period of 1979–85, driven by daily forcing of the same period. The close-to-a-century integration is certainly not long enough to guarantee any model to reach an equilibrium state. However, it is considered long enough to allow the models to demonstrate a clear pattern of how they evolve. Analyzing the integration results elucidates how the simulated Arctic ice–ocean system behaves in a relatively long integration, depending on whether and how it is constrained to the climatology, and yields useful information about polar climate modeling.
2. Model description
The model, with a horizontal resolution of 40 km and a vertical resolution of 21 levels of different thicknesses (see Table 2), is based on the ice–ocean model of Hibler and Bryan (1987). The ice model is of full viscous–plastic rheology (Hibler 1979, 1980) and is integrated using an efficient ice dynamics model that leads to a more rapid plastic solution for ice velocity (Zhang and Hibler 1997, hereafter ZH). The ocean model, based on that of Bryan (1969) and Cox (1984), is embedded with a mixed layer model of Kraus and Turner (1967) and has open boundaries at the Bering and Denmark straits and the Faroe–Shetland Passage for improving ocean circulation. Time-varying precipitation, together with a snow model, has been introduced to maintain a balanced salt budget and ocean stratification and to improve the system’s thermodynamics.
The basic setting for the coupled model is as follows. First, the sea–ice model is driven by atmospheric forcing, consisting of geostrophic winds, surface air temperature, humidity, and parameterized longwave and shortwave radiative fluxes. The ice model then supplies surface heat, salt, and momentum fluxes into the ocean as ocean surface boundary conditions. The ocean model, in turn, with the help of the mixed layer model, supplies current and heat exchange information to the ice model. The mixed layer model, driven by the surface fluxes, monitors the evolution of the mixed layer and modifies the temperature and salinity structure of the upper ocean. Some details of the mixed layer model and the embedding procedure are given in appendix A, and the ice model and the ocean model are briefly described below.
a. Sea–ice model
The sea–ice model is described in detail by Hibler (1979) and is only summarized here. The ice motion is determined by the following 2D momentum equation:
where m is the ice mass per unit area, u is the ice velocity, k is a unit vector normal to the sea surface, f is the Coriolis parameter, τa and τw are the air and water stresses, g is the acceleration due to gravity, ps is the surface pressure, and Fi is the ice interaction force due to internal ice stress. The air stress is determined by geostrophic winds, and the water stress is determined by the relative velocity between the ocean and the ice cover. The ice interaction force Fi is determined by using the nonlinear viscous–plastic constitutive law with an elliptic yield curve (Hibler 1979). The momentum equation is solved using the ice dynamics model of ZH, which improves the behavior of the plastic solution. The finite-difference scheme used with the momentum equation follows that of Hibler (1979), which ensures that the energy resulting from the solution of the momentum equation is bounded.
For mean ice thickness h and compactness A, the following continuity equations are used:
where and Γh and ΓA are thermodynamic source terms and Dh and DA are diffusion terms consisting of harmonic and biharmonic terms, which have been added to make the above ice advection equations numerically stable (Hibler 1979). Ice growth rate is calculated using the procedures of Semtner (1976) and Manabe et al. (1979), and the surface heat budget using computations similar to those of Parkinson and Washington (1979), as described in detail by Hibler (1980). Note that the sea–ice model consists of two categories, thin ice/open water and thick ice. In calculating the growth rate of the thick ice, a seven-level ice growth rate calculation (Walsh et al. 1985) is used. For each level of the thick ice, a weighted sum of the snow and ice conductivities is used (Riedlinger and Preller 1991).
The original ice model of Hibler and Bryan (1987) does not include a snow layer; it does, however, approximate the insulating effect of snow by allowing the ice surface albedo to be that of snow whenever the surface temperature is below freezing. We introduce a snow layer with a simple approach using monthly varying and spatially changing precipitation over the Arctic region (Vowinckel and Orvig 1970; Ranelli 1991). If the local ice condition is freezing, the precipitation is treated as snow accumulation; if it is melting or the precipitation encounters open water, the precipitation is treated as freshwater going directly into the ocean. The freshwater addition to the ocean is treated as a direct negative surface salt flux, whereas precipitation on the freezing ice, taken as snow, does not provide such a negative salt flux until it is melted at some point later.
b. Ocean model
where U is the horizontal ocean velocity, w is the vertical ocean velocity, ρ0 = ρw is the reference water density, and p is the pressure. The last term in (4) represents the surface momentum flux into the ocean calculated from the vector combination of air stress and ice interaction force (Hibler and Bryan 1987). The vertical and horizontal eddy viscosities are KM and AM, respectively;KM is chosen to be 500 cm2 s−1 in the top two levels and 1 cm2 s−1 elsewhere, and AM is uniquely chosen to be 8 × 107 cm2 s−1.
The temperature and salinity equations are
where T and S are, respectively, the ocean temperature and salinity, AH = 8 × 105 cm2 s−1 is the horizontal diffusion coefficient, and w′T′ and w′S′ are the vertical turbulent fluxes of heat and salt. Within the mixed layer, these flux terms are determined with the help of the mixed layer model (see appendix A for details). For the deeper layers below the pycnocline where the turbulence level is low, they are parameterized conventionally by eddy diffusivity such that
where KH is the vertical diffusivity coefficient chosen to be 5 cm2 s−1 in the top two levels and 0.2 cm2 s−1 elsewhere.
In (5) and (6), T0 and S0 are the climatological temperature and salinity taken from Levitus (1982), and Rt and Rs are the restoring constants. In the standard prognostic ocean simulation, Rt and Rs are taken as zero. In the first sensitivity study discussed here (denoted as ML-S hereafter), Rt is set to zero and Rs is set to 1/(30 days) within the mixed layer. In the second sensitivity study (denoted as Dp-S), Rt is set to zero and Rs is set to 1/(5 yr) under the mixed layer. In the third sensitivity study (denoted as Dp-ST), both Rt and Rs are set to 1/(5 yr). The first two sensitivity studies use diagnostic ocean models with different salinity-restoring terms. The basic concept of utilizing these models is to assimilate climatological ocean salinity data while still simulating seasonal and shorter-term salinity variations. At the same time, these two models allow the ocean temperature to vary freely without any constraint so that they can capture both long- and short-term temperature variations. However, the third sensitivity study allows only seasonal and shorter-term variations in both ocean salinity and temperature. Note that all except Dp-S correspond with some of the models listed in Table 1. Although Dp-S has not been used before, we feel that it has merit in demonstrating restoring effects on the ocean.
c. Numerical framework and forcing fields
The configuration of the model’s finite-difference grid and the bottom topography are shown in Fig. 1. The model has a horizontal resolution of 40 km × 40 km with a grid size of 130 × 102 in rectangular coordinates. The grid is divided into three geographic areas, as shown in the figure, for the purpose of analysis. The ocean’s vertical dimension has 21 levels of different thicknesses, as shown in Table 2. Bottom topography, from the ETOPO5N database at the National Geophysical Data Center (see Ranelli 1991 for details), is resolved by using differing numbers of levels at different locations. For very shallow areas, four levels are specified, which means the minimum ocean depth is 79.4 m.
Besides using nonslip lateral boundary conditions along the ocean walls, the model uses open boundary conditions across the Bering Strait, Denmark Strait, and Faroe–Shetland Passage. A water inflow of 0.8 Sv (1 Sv ≡ 106 m3 s−1) is specified at Bering Strait, an inflow of 5.6 Sv at most of the Faroe–Shetland Passage, and an outflow of about 6.4 Sv at Denmark Strait and part of the Faroe–Shetland Passage (see Table B1). In addition to the specified mean flow across the open boundaries, a geostrophic adjustment is used to achieve vertical velocity variations due to local density variations. Appendix B gives a detailed description of the treatment of the open boundaries.
As mentioned above, all the models used in this study were run repeatedly for ten cycles over the 7-yr period of 1979–85. In the first seven cycles, a 6-h time step was used in the integration of the ocean conservation equations, whereas a 0.3-h time step was used for the ocean momentum equations, a so-called distorted physics (Bryan 1984). The one-quarter-day time step was also used for solving equations in the ice and mixed layer models. In the last three cycles, synchronous simulations were carried out to make the solution dynamically consistent. The time step for the synchronous integration was set to 0.3 h for all equations. Solution was based on a leapfrog time-stepping procedure, but a modified Euler time step was used after every 17 leapfrog time steps to eliminate one branch of the split solution caused by leapfrog differencing.
The model is driven by daily geostrophic winds, which were derived from daily sea level pressure fields obtained by merging the analysis of Arctic buoy data with the NCAR analysis of daily sea level pressure data (J. E. Walsh 1991, personal communication). In addition, Hansen’s dataset of monthly mean air temperature (Hansen and Lebedeff 1987) is used with some modifications to take into account the effect of ice feedback on air temperature (see Hibler and Zhang 1993). The shortwave and longwave radiation is calculated using the formulation of Parkinson and Washington (1979). The river runoff is the same as that specified by Hibler and Bryan (1987), and the monthly varying, region-dependent precipitation is from Vowinckel and Orvig (1970) (see Table 3). Shown in Fig. 2 are the 7-yr (1979–85) mean downward longwave and shortwave radiative fluxes, the surface air temperature, and the geostrophic winds, which are used in this study and a number of previous studies (Ranelli 1991; Hibler and Zhang 1993, 1994; Zhang 1993; Ip 1993; Flato and Hibler 1995).
3. Results and discussion
For clarity, we list below the standard prognostic run (denoted as Std) and the major sensitivity runs
Std: no restoring, standard prognostic run
ML-S: 30-day climate restoring of the mixed layer (ML) salinity
Dp-S: 5-yr climate restoring of the salinity below the ML
Dp-ST: 5-yr climate restoring of the salinity and temperature below the ML.
We also conducted a fourth sensitivity study, ML-ST, which has 30-day climate restoring of both the mixed layer salinity and the mixed layer temperature. However, this model overmelts the ice because constraining the temperature to the warm Levitus values generates excessive heat in the mixed layer. Therefore, its results are not presented here, but some comments about this model are given in section 4. As mentioned earlier, the standard run and all the sensitivity studies were first integrated asynchronously for seven cycles over the 7-yr period of 1979–85. Then three additional cycles of synchronous integrations were conducted. The initial oceanic state was set at the Levitus (1982) climatological temperature and salinity. The initial ice thickness and concentration fields were obtained after a preliminary run of 7 yr with the standard model. Most of the results shown below are based on the tenth 7-yr simulation period unless stated otherwise.
a. Ocean salinity
In this section and some following sections, the ocean is divided into two layers: the “Upper Layer,” defined as the top four ocean levels which goes as deep as about 80 m, and the “Lower Layer,” defined as the rest of the ocean. This was done to simplify analysis of the oceanic processes, given that the simulated depth of the mixed layer is generally close to the 30–50 m estimated from observations (Coachman and Aagaard 1974) and is contained within the top four levels of the model. Figure 3 shows the evolution of the simulated annual average salinities of the Upper Layer and Lower Layer for the Arctic Ocean in the course of the first 49-yr asynchronous integrations and then the 21-yr synchronous integrations. Also shown are the corresponding Levitus values. As can be seen from the figure, after a few 7-yr cycles of asynchronous simulations, the salinities of both the Upper and the Lower Layers are all relatively stabilized, although they may not reach equilibrium. When the simulations are changed to be synchronous, the salinities are changed accordingly and start to move toward steady state again. This indicates that these models are well behaved in terms of adjusting themselves to maintain a salt balance.
In the Upper Layer of the Arctic, all the models except ML-S predict an increase in salinity compared with the Levitus value. This is because, for Std, Dp-S, and Dp-ST, which do not have salinity constraints in the mixed layer, there is a large amount of salt rejection from the ice cover. As a result, even when the contributions of negative salt flux due to precipitation and river runoff and the effects of other internal oceanic processes are taken into account, the ocean has a net gain of salt in the Upper Layer until reaching a relatively steady state. The lower salinity predicted by ML-S is evidence of the constraints to climatology in the mixed layer even though a large amount of salt rejection still occurs at the ice–ocean interface. Note that ML-S, unlike Dp-S and Dp-ST, does not prevent the state of the deeper Arctic water from shifting far away from the climatology (climate drift). ML-S does, however, exercise more control in the mixed layer, via the mixed layer salinity corrections, so that the mean salinity predicted in the mixed layer is closer to the climatological mean as shown later. But why would ML-S predict a much lower mean salinity than the Levitus value in the upper layer? The reason is that the incoming salt at the surface due to ice formation is immediately taken away by the restoring term incorporated in the mixed layer and hardly reaches deeper layers. As a result, the water between the base of the mixed layer and 80-m depth and in even deeper layers is fresher than the Levitus value (see Fig. 4). In fact, all models predict lower salinity in the lower layer than the Levitus values (also see Fig. 4).
However, Dp-S and Dp-ST are rather effective in maintaining a state that is closer to the initial climatological state because of their 5-yr weak salinity restoring in the layers below the mixed layer. Introducing ocean temperature restoring in Dp-ST does not cause a significant difference in salinity. Therefore, as diagnostic models, Dp-S and Dp-ST are more suitable for a short-term study as far as salinity is concerned. The other two models predict more deviations from the initial condition and take a longer time to reach a relatively steady state. Again, the salinity restoring in the mixed layer affects the salinity in the layers below so that the Lower Layer salinity predicted by ML-S is considerably lower than that predicted by the other models; in contrast, the free salinity evolution of Std results in a salinity that is moderately lower than the Levitus value, which is also illustrated in Fig. 4.
Figure 4 shows distributions of Levitus salinity and salinity differences between model predictions and Levitus data on a transect that crosses the Arctic Basin and connects Alaska and Franz Josef Land (see Fig. 1). As can be seen from Fig. 4a, Std predicts a higher salinity than the Levitus data in the upper 200 m in much of the Arctic except near Franz Josef Land, where it predicts a lower salinity. Below 200-m depth, there is no significant drift from the Levitus data. The difference between the salinity simulated by ML-S and the Levitus data is shown in Fig. 4b. Clearly, ML-S does a good job of maintaining a surface salinity close to the Levitus values, but right below the surface layer, it predicts water that is considerably fresher than the Levitus values. This means that the 30-day salinity restoring terms used in the mixed layer are rather effective in taking away the salt rejected from growing ice. This is why the mean salinities in both the Upper Layer and the lower layer are much lower than the Levitus values (Fig. 3). (The salinity distribution of Dp-ST is very close to that of Dp-S; we therefore show only Dp-S results in Fig. 4.) In the layers below about 100-m depth, the salinity predicted by Dp-S is close to the Levitus data, indicating that the 5-yr salinity-restoring terms are good enough to constrain the salinity to climatology in the deeper layers.
The features shown in Fig. 4 are further illustrated in Fig. 5, which includes the mixed-layer salinity fields simulated by Std, ML-S, and Dp-S and the Levitus surface salinity field for comparison. Again, the field derived by Dp-ST is similar to that derived by Dp-S and is therefore not shown. There are major differences in the three fields predicted in the Arctic. The surface water derived by Std is more saline than that of the Levitus data owing to salt rejection from the ice cover. The surface water derived by ML-S is close to the Levitus values but is probably excessively fresh since most of the Levitus data were collected in the summer. The surface salinity derived by Dp-S, and Dp-ST as well, is somewhat in between. This is because the 5-yr constraint on the layers below the mixed layer allows the surface salinity to be more influenced by the surface salt flux, which results in a considerable change in the surface salinity compared with the initial conditions set by the climatology.
Why would the prognostic model, Std, predict considerably saline surface water in the central Arctic? Because Std does not have the restoring term for salt correction, its mixed-layer salinity field is closely related to the surface salt flux and hence to the ice growth rate shown in panel 1 of Fig. 6. Also included in Fig. 6 are the ice growth fields simulated by the other three models, which will be discussed later in section 3d. Although for most of the central Arctic, the predicted annual ice production is not far away from the 0.9 m estimated by Maykut (1982), the ice growth tends to dump considerable amounts of salt into the ocean annually, as shown in the left panel of Fig. 7. The salt, in turn, tends to be trapped in the gyre of the Beaufort Sea and Canadian Basin, and after sufficiently long prognostic integration, the model eventually predicts a high salt concentration in the central Arctic. This high salinity in the central Arctic has consequences for the ocean surface circulation, as mentioned in section 3c(1). The right panel of Fig. 7 shows the surface salt flux predicted by ML-S due to both ice growing/melting and the salinity correction terms in the mixed layer. Note that the surface salt flux due to the salinity correction is large enough to overwhelm that due to ice growing/melting. As can be seen from Fig. 7, there is a large negative salt flux in the Canadian Basin and the eastern Siberian Sea, and a large positive salt flux in the Eurasian Basin and Laptev Sea. How would this pattern be developed? The negative salt flux in the Canadian Basin may be attributed to more salt rejection from the ice cover in that region (see upper right panel of Fig. 6). The ocean surface circulation pattern may also play a role in bringing fresher mixed layer water from the Canadian Basin to the Eurasian Basin, which triggers the restoring terms to generate a positive salt flux in Eurasian Basin.
b. Ocean temperature
Figure 8 shows the time series of the simulated average ocean temperature in the Upper and Lower Layers of the Arctic Ocean along with the Levitus value. There is a clear shift from the asynchronous mode to the synchronous mode in the Lower Layer results but not in the Upper Layer results. The model predictions of the average temperature in the Upper Layer of the Arctic show little difference; they are all much smaller than the corresponding Levitus value, which probably indicates, again, that the climatology has a certain bias in favor of summer conditions. The slightly higher temperature predicted by Dp-ST is due to the temperature correction imposed below the mixed layer.
In the Lower Layer, however, there are some major differences. The temperature predicted by ML-S is much higher than that predicted by the other models and is drifting rapidly away from climatology, whereas that predicted by Dp-S is decreasing continuously. This indicates that the heat balance in the Arctic Ocean’s deeper layers is rather sensitive to whether and how salt-correction terms are introduced in the models. As pointed out later, climate restoring has a great impact on both vertical and horizontal oceanic heat transfer once the ocean temperature is allowed to be free, since without an interactive atmosphere on top, an ice–ocean system is not a closed system. In this open system, heat gained or lost to the atmosphere is not reversible. Therefore, any excessive gain or loss of heat in the ocean may result in a sustained temperature increase or decrease. However, this does not necessarily mean that ML-S or Dp-S would never reach a thermal equilibrium; rather, it tells that a long integration is needed for them to reach such an equilibrium, which will be far away from the climatological state.
Although it is impossible for Std to reach an equilibrium in a 70-yr simulation, the mean temperature predicted in the Lower Layer does look roughly like it is in a steady state and is not far from the climatology. This may indicate that without the interference of the climate restoring terms, the standard prognostic model is able to adjust the related oceanic processes and to maintain a rough heat balance in the 70-yr integration period. We do not know, however, what would happen if this model were integrated for thousands of years, an integration beyond the capacity of our computers. Therefore, we are content with the 70-yr simulation results, which at least show the trends in the models’ drift relative to the initial climatological state. The results show that the thermal state predicted by ML-S and Dp-S tends to drift away from the climatology at a much faster rate than the drift, if any, predicted by Std. As for Dp-ST, its mean temperature quickly reaches a steady state after about one 7-yr cycle and stays the closest to the climatology. Therefore, Dp-ST is more suitable for short-term integration if one is interested in short-term variations in the ocean density structure and ocean circulation from some climatology (e.g., Hibler and Zhang 1994). Its predictability of long-term variations, however, is limited.
Why does the Lower Layer temperature predicted by ML-S increase more rapidly than that predicted by Std in the very beginning of the simulation and keep climbing over the whole 70-yr integration period? Part of the answer can be found in Fig. 9, in which predictions of the vertical heat loss in the Arctic Lower Layer due to different oceanic processes are plotted for the early stage of the simulation (cycle 1) and the late stage of the simulation (cycle 10). ML-S predicts much less convective overturning in cycle 1 because the salt flux correction in the mixed layer changes the buoyance forcing at the surface; as a result, the heat loss due to this process is considerably reduced compared to that predicted by Std. Meanwhile, the heat loss due to vertical advection is also considerably smaller than that predicted by Std because the divergence of the upper-ocean velocity is different. This means the upper-ocean dynamics is different and it is found in section 3c to be closely related to horizontal density distribution and vertical stratification. Because of the inability of these two oceanic processes to deliver enough heat upward, ML-S traps more heat in the intermediate layers of the ocean, so the temperature of the Lower Layer tends to increase from the very beginning. Once the intermediate layers become overly warm, after a sufficiently long integration such as that represented in cycle 10, the limited overturning process in ML-S begins to deliver more heat upward, compared with Std, in an attempt to make the heat balance. The attempt is quite successful in that it brings the total vertical heat loss predicted by ML-S close to that predicted by Std.
Because the total vertical heat loss predicted by Std and ML-S in the last cycle is about the same, why does the ocean temperature formed by ML-S keep climbing while that of Std stays relatively flat? What happens is that ML-S traps heat not only in the Arctic Basin, but also in the ice-covered area of the Greenland Sea and, particularly, in the Fram Strait area. This is because these areas, like the Arctic, are also deep and ice covered and thus more subject to the influence by the surface salinity restoring. In these areas, the Lower Layer temperatures predicted by all models except Dp-ST have the same behavior shown in Fig. 8. Figure 10 shows the temperature distributions predicted across Fram Strait by the three simulations Std, ML-S, and Dp-S without temperature corrections. The temperatures predicted by ML-S are much higher than those predicted by the other two models as well as those of the Levitus data. The Std results seem reasonable in that they reflect a layer of warm Atlantic water coming into the Arctic as the West Spitzbergen Current. The Dp-S temperatures are the lowest of the three cases, while the warm Levitus temperature data are expected because of their summer bias. Since ML-S simulates an Atlantic water inflow of about 2.0 Sv at Fram Strait [see Fig. 16 in section 3c(3)], we expect that more heat is carried into the Arctic than predicted by the other two models. This is shown in Fig. 11, which compares model predictions for the heat budget in the Arctic Ocean below 80 m in the tenth cycle. There are major differences in the heat inflow at Fram Strait because of the different temperature distributions and strengths of the Atlantic water inflow (Fig. 16); in contrast, the differences in vertical heat loss are small among the three models without temperature corrections (Std, ML-S, and Dp-S). ML-S predicts much more heat inflow at Fram Strait and hence an excessive net heat gain. This is why the temperature in the ocean interior keeps climbing even though the vertical heat loss is close to that simulated by Std. In contrast, Dp-S predicts a net heat deficit owing to the low temperature and low Atlantic water inflow at Fram Strait; Std predicts a smaller heat gain than ML-S and Dp-ST, which, although not in an exact heat balance, is nevertheless close. All this discloses that, for the models with freely evolving temperature, the lateral heat transport is crucial in determining the heat budget during the later stage of simulation. It is, however, the vertical processes in the course of spinup that trigger such an outcome by trapping heat in the intermediate layers. The role that the vertical processes play in the later stage also should not be neglected. The overturning process in ML-S, for example, shows a tendency to adjust itself and release more heat upward in an attempt to create a heat balance.
Because of the temperature constraints, Dp-ST behaves quite differently from the other three models in terms of heat budget, as shown in Fig. 11. The Lower Layer of Dp-ST gains substantially more heat vertically in the summer and loses more heat in the fall. Because of the excessive summer heat gain, the Lower Layer has a net heat gain of about 7 × 1012 W. This heat surplus is obviously “corrected” by the temperature restoring below the mixed layer. This is why the Lower Layer mean temperature predicted by Dp-ST is in a steady state as shown in Fig. 8. Why would the Lower Layer of Dp-ST gain more heat in summer and lose more in fall? This is because the related oceanic processes have peculiar behaviors, which are shown in Fig. 12. This figure is a replot, on a different scale, of the first three panels of the right column of Fig. 9 with the Dp-ST results added. Figure 12 shows that the Lower Layer of Dp-ST loses more heat in fall by convective overturning. The reason is that in fall the ocean surface layer starts to cool owing to atmospheric forcing, while the deeper layers are kept relatively warm by the temperature-restoring terms, which induces more overturning that brings up more heat. In summer, in contrast, the surface water is warm, while the deeper layers are kept relatively cool by the restoring terms; therefore, the vertical diffusion process delivers considerably more heat to the deeper layers. This is why the lower left panel of Fig. 11 shows that in the Lower Layer Dp-ST predicts considerable vertical heat gain, which is “corrected” by the temperature restoring.
c. Ocean velocity
1) Ocean surface velocity
Figure 13 shows simulated 7-yr mean fields of ocean surface velocity, defined as those at the third level of the ocean (about 37 m deep), which roughly represent the surface geostrophic velocity. Std predicts a clockwise gyre over the Beaufort Sea and Canadian Basin and a transpolar drift stream that is probably a little too weak. There is also a strong flow along the shelf break off Franz Josef Land and Spitzbergen. Outside the Arctic Basin, the model captures some major surface currents, such as the East Greenland Current, the Norwegian Current, the North Cape Current, the East and West Spitzbergen Currents, and the East Icelandic Current. ML-S predicts much stronger surface circulation, with a clockwise gyre that covers almost the whole Arctic Basin. The transpolar stream moves more toward the Canadian Basin instead of the Greenland Sea, and the flow along the shelf break off Franz Josef Land and Spitzbergen is much stronger. Dp-S and Dp-ST predict a stronger transpolar stream but a weaker flow along the shelf break off Franz Josef Land. The circulation gyre in the Arctic tends to be limited to the Beaufort Sea area, and the flow in the Greenland Sea is weaker than predicted by the other two models. Between Dp-S and Dp-ST, there is no significant difference in the circulation pattern in most of the Arctic, whereas there are some differences in the Greenland Sea. This indicates that the temperature restoring below the mixed layer does not significantly affect the upper circulation except in the Greenland Sea, where the coming warm Atlantic water is more likely to be modified by the restoring terms.
It seems clear that all the differences in the upper-ocean circulation predicted by these four models are more closely related to the spatial salinity distribution in the upper ocean than to the temperature distribution;this, in turn, is affected by the surface salt flux conditions and salinity restoring schemes as mentioned earlier. With the dome-shaped surface salinity distribution in the central Arctic, predicted by Std, the baroclinic component of the flow tends to form an anticlockwise circulation in the upper layer, which is in the direction opposite to that induced by wind curl or ocean surface stress. On the other hand, when the salinity in the central Arctic is lower than that in the surrounding area, then the flow component driven by buoyancy and the component driven by wind curl or ocean surface stress are in the same direction, and a strong clockwise gyre is created. This means that buoyancy plays an important role in driving the circulation, and this is why ML-S predicts the strongest surface flow in the Arctic, followed by Dp-S, and Std, while the strength of the surface flow predicted by Dp-ST is close to that predicted by Dp-S. The magnitude of the ocean surface velocity predicted by Std may be a little underestimated owing to Std’s higher surface salinity, whereas that predicted by ML-S may be overestimated. On the other hand, Dp-S and Dp-ST, because of their intermediate surface-water salinity, may be more realistic in simulating a moderate surface flow and a stronger transpolar stream.
2) Ocean velocity at an intermediate layer
To show the behavior of the ocean velocity at a deeper layer, Fig. 14 illustrates the velocity fields in level 9, at a depth of about 400 m, which is close to the core of the Atlantic water inflow. The panels all show a velocity gyre in the Canadian Basin and Beaufort Sea, but the gyres simulated by Std and ML-S are much stronger. Elsewhere in the Arctic, Std predicts a narrow anticlockwise gyre over part of the Eurasian basin, whereas Dp-S and Dp-ST predict a gyre that covers a larger area. All predict strong Atlantic water inflow across Fram Strait near the west coast of Spitzbergen. The circulation pattern predicted by Dp-S and Dp-ST is more influenced by the bathymetry and therefore is probably more realistic. In particular, the Dp-ST results clearly show the flow steering effects of the Lomonosov, Alpha, and Northwind Ridges, the Arlis Plateau, and the Chukchi Cap. At Fram Strait, the Dp-ST results show a more distinctive Atlantic water inflow in the West Spitzbergen Current and a cold Arctic water outflow in the East Greenland Current compared with Dp-S results. This indicates that the temperature restoring below the mixed layer does influence the flow pattern in the intermediate layers and probably makes the flow pattern look more realistic. Also note the differences in the strengths and patterns of the circulation predicted by the models for the GIN Sea.
3) Flows across Fram Strait
The velocity field generated by Std in Fig. 13 and a calculation of integrated water-transport streamfunction indicate that the 5.6 Sv of Atlantic water specified in the simulation flows north into the GIN Sea through the Faroe–Shetland Passage. Most of the water (about 3.3 Sv) seems to turn around to join the East Greenland Current and flow back into the Atlantic. About 0.8 Sv flows farther northward into the Barents Sea, while about 1.2 Sv tries to force its way into the Arctic Basin through Fram Strait. The latter is not very successful at the surface level because it encounters the massive cold Arctic surface water flowing out of Fram Strait. Therefore, this water often dives to an intermediate level, where it flows farther into the Arctic Ocean. This scenario is supported by Figs. 15 and 16.
Figure 15 shows the velocity distributions on a transect across Fram Strait. Again, as in Fig. 14, Atlantic water flows into the Arctic mainly at the east side of the strait. The depth of the inflow ranges from the surface to about 1000 m deep, with the core being centered at about 400 m. In contrast, the massive outflow of cold, fresh Arctic water mainly occurs above 400-m depth, and the saline Arctic water mainly flows out below 1000-m depth. ML-S predicts more Atlantic water inflow at the west side of Fram Strait, which pushes the outflowing cold Arctic water more toward the surface there. Dp-S predicts the smallest inflow of Atlantic water as well as the smallest outflow of cold, fresh Arctic water, which is probably less realistic. However, Dp-S predicts more outflow of Arctic saline water below 1000-m depth. The Dp-ST results are interesting in that they show the strongest Atlantic water inflow at the east side of Fram Strait, with the depth of the flow ranging from the surface to as deep as 1800 m. On the west side of Fram Strait, there is no Atlantic water inflow at all, only the Arctic water outflow. This kind of flow pattern at Fram Strait is probably more realistic. Again, one can see that the deep temperature restoring has more impact in the Greenland Sea and Fram Strait area because of the large lateral heat advection there.
Figure 16 shows the model predictions for the annual mean inflow and outflow at Fram Strait. Note that the volume of Atlantic water inflow predicted by the models is close to the 0.9 Sv estimated by Rudels (1987). The standard model puts the inflow of Atlantic water in the neighborhood of 1.2 Sv and the outflow of Arctic water in the neighborhood of 2.8 Sv, with some interannual variations. Both ML-S and Dp-ST predict about 0.7 Sv more inflow and 0.5–0.8 more outflow than Std; Dp-S predicts slightly less inflow and slightly more outflow. This indicates that the climate restoring affects the communication between the Arctic and the Greenland Sea in a two-way manner. Inside the Arctic Ocean, climate restoring changes the density structure, which affects the strength of the water mass transport and hence the outflow of Arctic water. The more Arctic water that flows out, the more Atlantic water that is likely to flow in, and vice versa. On the other hand, the restoring would also change density structure of the GIN Sea, thus changing the amount of Atlantic water inflow and therefore the amount of Arctic water outflow. For the deeper outflow of the saline Arctic water, Dp-S predicts a rate of 1.5 Sv, which falls into the range of 1.3–2.0 Sv estimated by Aagaard et al. (1991); Dp-ST predicts a rate of 0.88 Sv, which falls into the range of 0.7–1.0 Sv estimated by Heinze et al. (1990), while both Std and ML-S predict an outflow of only 0.54 Sv, which may be a little low.
Among the four models, which is the “best” in terms of predicting ocean circulation? This is of course a tricky question. It seems clear, however, that ML-S gives the poorest results. The surface currents from this simulation exhibit a transpolar drift stream that is amplified not over the North Pole, but, in fact, near the continental shelf break off Franz Josef Land and Spitzbergen. Also, the currents in the intermediate layers predicted by ML-S are the most decoupled from bathymetry of all the models, in contrast to recent icebreaker and submarine observations (Anderson et al. 1994; Morison et al. 1998). In Fram Strait, the deep outflow is weaker than observed, and the less-than-realistic intermediate inflow in the East Greenland Current predicted by all models (except Dp-ST) is strongest in ML-S. The failure of the model indicates simply that the surface (Levitus) salinity field has errors, which is to be expected given the poor data coverage in the Arctic Ocean (Levitus 1982). The situation is not appreciably improved in the newer version of this dataset (Levitus 1994).
Although the results of Std are clearly better than those of ML-S, the problem is that Std predicts a high salinity dome in the Canadian Basin that weakens the surface currents so that the surface transpolar stream is rather weak. The surface circulation in the Beaufort Sea may be weak too. In the intermediate layers, the flow is also more or less decoupled from the bathymetry. At Fram Strait, the Atlantic water inflow is more toward the east side of the strait, which is relatively realistic compared with the results of ML-S. The deeper outflow, however, is also weak, as it is with ML-S. Still, it is encouraging that the ocean currents simulated by Std are reasonably close to reality, given that the model does not have climatological corrections.
Probably the best ocean circulation pattern is generated by Dp-ST. Since both salinity and temperature are constrained to climatology, the results are more like those “observed,” which is certainly not surprising. Because DP-ST predicts a relatively flat salinity field in the central Arctic, the surface circulation pattern has a nicely shaped transpolar stream and Beaufort gyre, which matches reasonably well the ice drift pattern observed in most of the Arctic (Rigor 1992; Pfirman et al. 1997). Furthermore, the flow pattern of the intermediate water clearly shows bathymetric steering, and the Atlantic water inflow is concentrated mostly at the east side of Fram Strait, which is probably most realistic. Although Dp-S generates about the same surface circulation in the Arctic as Dp-ST, its prediction of intermediate water flow is clearly inferior to that of Dp-ST. The Atlantic water inflow predicted by Dp-S at Fram Strait is probably too weak.
d. Influence on ice
Can these different models, which directly affect the simulation of oceanic salinity and temperature, have an influence on the prediction of sea ice? The answer is yes because, as shown before, they affect the ocean dynamics and thermodynamics, which, in turn, affect the ice conditions. In this section, we examine how the ice cover is influenced. To obtain an idea of how the ice dynamics is affected, we plotted the ice motion field predicted by Std and the difference between it and the fields predicted by ML-S, Dp-S, and Dp-ST. The results are shown in Fig. 17. In the Arctic, the ice motion pattern predicted by Std (Fig. 17a) consists of a circulation gyre and a transpolar stream. This pattern is in an agreement with the observed drift of Arctic buoys (Rigor 1992; Pfirman et al. 1997). ML-S (Fig. 17b) predicts considerably more ice motion than Std owing to stronger surface currents, resulting in a circulation gyre that covers the entire Arctic Basin. Although there are some differences between the ice motion predicted by Dp-S and Std and between Dp-ST and Std, they are smaller than the difference between ML-S and Std because the difference in ocean surface velocity is relatively smaller (see Fig. 13). The transpolar motion and the Beaufort gyre predicted by Dp-S and Dp-ST are larger than those predicted by Std, whereas the ice drift along the shelf off Franz Josef Land and Spitzbergen is smaller, owing again to the surface ocean circulation. To evaluate the overall behavior of the simulated ice drift, we conducted a statistical analysis (Zhang 1993) based on the ice velocities simulated by the four different models and all the buoy-drift data available in the whole region over the 7-yr period of 1979–85. The buoy data show a mean drift of 6.24 km day−1; the simulated drifts are 6.83 km day−1 (Std), 7.43 km day−1 (ML-S), 6.77 km day−1 (Dp-S), and 6.87 km day−1 (DP-ST). Thus statistically, the ice simulated by Std drifts about 9% faster than the observed drift, that by ML-S about 19% faster, that by Dp-S 8% faster, and that by Dp-ST 10% faster. Clearly, ML-S behaves the worst in computing ice velocity, while the other models do not differ significantly. Note that the air stress on the ice is calculated by bulk parameterization, with the stress coefficient and turning angle being estimated based on summer observations (McPhee 1980), which tends to slightly overpredict ice motion.
As a further comparison, we plotted (Fig. 18) the simulated seasonal ice concentrations in the Arctic and the observed ice concentrations. The observed results are from the Satellite Multichannel Microwave Radiometer (SMMR), and the simulated results are from the 10th cycle of model integration. The models often slightly overpredict ice coverage in winter and underestimate it in summer. The winter ice concentration predictions are rather insensitive to the model used, whereas the summer concentrations exhibit some differences among the models. Overall, the summer concentrations predicted by Std and Dp-S are the closest to the observations, whereas ML-S and Dp-ST tend to underpredict the ice concentrations. The Std and Dp-S results are close to each other, which is conceivable because both models are free of constraints on the thermal state. However, the ML-S and Dp-ST results also seem to stay close to each other. The reason is not clear. One possibility could be that the Dp-ST’s temperature correction terms below the mixed layer play a role similar to that of the overwarm intermediate layers of ML-S in reducing the ice extent.
To obtain an idea of how the spatial ice thickness distribution is influenced by the models, we plot the ice thickness field predicted by Std and the differences between the thickness fields predicted by ML-S and Std, by Dp-S and Std, and by Dp-ST and Std (Fig. 19). The pattern of the thickness field predicted by Std (Fig. 19a) is in reasonably good agreement with the pattern observed by LeShack (shown by Hibler 1980), with the predicted field having a slightly smaller ice thickness. Compared with Std, ML-S and Dp-ST predict a 10%–25% thinner ice cover in most of the Arctic (Figs. 19b and 19c). The reason ML-S predicts thinner ice is that it generates a stronger circulation at the ocean surface and thus tends to move more ice out of the central Arctic toward Franz Josef Land and Spitzbergen. Meanwhile, there is substantially more oceanic heat flux in the Eurasian basin and in the Franz Josef Land and Spitzbergen area, as shown in Fig. 20b, which melts more of the ice growing in or coming into that area. This scenario is also reflected in the upper right panel of Fig. 6, which shows considerably less ice growth in those areas. This is why ML-S predicts less ice in the central Arctic, the Eurasian basin, and the Franz Josef Land and Spitzbergen area.
Why would ML-S generate such a large oceanic heat flux in the Eurasian basin and in the Franz Josef Land and Spitzbergen area? This is because in those areas ML-S creates a larger lateral heat advection owing to the excessive heat carried in by the incoming Atlantic water before it mixes with the Arctic water, as shown in the top panel of Fig. 11. Moreover, the intermediate layer formed by ML-S is especially warm so that more heat tends to be transported into the mixed layer from the warm intermediate layer by vertical processes such as convective overturning and vertical advection. This underlines the relevance of ocean dynamic and thermodynamic processes to the Arctic ice conditions. As pointed out before, the excessive heat carried in by the incoming Atlantic water is owing to a warm intermediate layer in the Fram Strait area and the warm layer is, in turn, owing to the suppressed vertical convective overturning and advection in the course of model spinup (Fig. 9). This further shows the importance of correctly simulating the vertical overturning and advection processes in both the Arctic and adjacent seas. If these two processes are not reasonably represented, the ice-covered ocean in the Arctic and around Fram Strait may become too warm, and the resulting excessive heat carried into the Arctic by the Atlantic water inflow may gradually shrink the ice sheet there over a long simulation time. Note that in the central Arctic and Beaufort Sea, ML-S creates slightly less oceanic heat flux than Std. This is because the upper layer predicted by ML-S is much fresher than that by Std, which reduces the level of convective overturning and upward heat transport even though ML-S has a warmer intermediate layer. This is why the basin-wide-averaged vertical heat loss at 80-m depth from ML-S and Std is close, as shown in the right bottom panel of Fig. 9.
The reason for Dp-ST predicting a shrunken ice pack is different from that for ML-S. Because Dp-ST incorporates temperature restoring terms, the variations in the model’s thermal state are constantly “corrected” by a three-dimensionally distributed heat source/sink below the mixed layer. This inevitably leads to more oceanic heat flux from the deep ocean due to convective overturning in fall, as shown in the top panel of Fig. 12. This is why we see more oceanic heat flux in the areas of the Beaufort and Chukchi Seas, the Eurasian basin, and Franz Josef Land and Spitzbergen in Fig. 20d. Particularly, there is substantially more ocean heat flux in the Fram Strait and Spitzbergen area, which again indicates that the Levitus temperature data there are particularly warm, as shown in the bottom panel of Fig. 10. Because of this, the ice growth rate predicted in the Arctic is a little smaller than the rates predicted by Std and Dp-S, as shown in the lower right panel of Fig. 6, and hence the Arctic ice pack is smaller. In summer, the vertical diffusion process generated by Dp-ST would transfer more heat to the deeper ocean layers from the surface, as shown in the bottom panel of Fig. 12. That part of the heat may not affect the ice growth and decay because it would be “corrected” by the restoring terms.
Compared with Std, Dp-S predicts only slightly less ice in the central Arctic but substantially more ice in the area north of Franz Josef Land and Spitzbergen, as shown in Fig. 19c. This is because Dp-S forms a slightly stronger transpolar stream of water and ice (see lower left panel of Fig. 13 and Fig. 17c), which tends to drive more ice toward Fram Strait so that less ice is left in the central Arctic and the Beaufort Sea and more ice enters the Spitzbergen area and the Greenland Sea. Unlike ML-S and Dp-ST, the ice entering the Spitzbergen area and the Greenland Sea does not encounter a large oceanic heat flux and is therefore not melted. This is why Fig. 19c shows more ice in the Greenland Sea.
The different ice thickness distributions predicted by these models are an accumulative effect of ocean circulation on the ice dynamics and thermodynamics. Another measure of the accumulative effect is seen by examining the ice outflow at Fram Strait (Fig. 21). The amount of ice moving out of the strait is of interest because obviously it is closely related to the salt/freshwater budget in the Arctic. On the other hand, the Arctic ice outflow may also affect the thermohaline circulation in the Atlantic and probably in the global ocean (see Hibler and Zhang 1995). Also plotted in Fig. 21 is the areal transport estimated by Colony (1990) based on buoy observations. Compared with the observational estimate, all the models somewhat overestimate the areal transport. The accumulative effect of ocean circulation on the ice drift is reflected in the slightly higher areal transport predicted by ML-S, which forms particularly robust surface currents in the Arctic, particularly in the Franz Josef Land and Spitzbergen area. However, its volume transport is lower than that of Std, as is that of Dp-ST, because of a thinner ice cover. In contrast, the volume transport of Dp-S is higher than that of Std because of a slightly stronger transpolar stream that pushes thicker ice out of the Arctic, as mentioned before.
4. Concluding remarks
This study has two objectives. One is to construct a fully prognostic ice–mixed layer ocean model of fine resolution for the Arctic, Barents, and GIN Seas and to examine the characteristics of the model in order to find a way to improve the numerical predictability of the coupled Arctic ice–ocean system without constraining it to its historic mean state. The other is to understand the implications of climate restoring, as commonly used by the Arctic modeling community, on the simulated Arctic ice–ocean circulation by comparing the results of the standard prognostic model with three major sensitivity studies using different diagnostic models. All the model integrations are based on a 70-yr period moving from asynchronous simulations to synchronous simulations at year 50. As pointed out before, the integration period is not long enough to guarantee an equilibrium oceanic state, but it is enough to allow the models to show a trend in how the modeled oceanic state and the sea ice pack evolve with and without climatological corrections. Although we do not know how each model would drift relative to climatology in an even longer integration such as thousands of years, we found that the close-to-a-century integrations yield considerable information about modeling the Arctic ice–ocean system with and without climate restoring.
The standard prognostic model basically captures the main features of the ice and ocean circulation in the Arctic and the adjacent areas. The ice motion, concentration, and thickness are in reasonably good agreement with observations. The model reproduces some major ocean currents at the surface that are close to those inferred from density distributions and long-term station drifts (see Fig. 17 of Coachman and Aagaard 1974). In the Arctic, particularly, the model captures the clockwise ocean circulation gyre in the Beaufort Sea and Canadian Basin area and the transpolar drift stream in the central Arctic and Fram Strait area. However, the Beaufort Sea/Canadian Basin gyre and the transpolar stream may be weaker than in reality owing to the unrealistically dome-shaped high salinity distribution in the central Arctic. Outside the Arctic, the modeled surface currents seem reasonable [see section 3c(1)]. In the intermediate layers, the model predicts a narrow and weak anticlockwise circulation in the Eurasian basin, as is also simulated by a diagnostic model used by Aukrust and Oberhuber (1995). However, the robust flow does not signal a strong influence by the bathymetry. The inflow of Atlantic water at Fram Strait occurs mainly west of Spitzbergen in the intermediate layers above 1000 m, which carries about 1.2 Sv of water into the Arctic. This value is close to the 1 Sv estimated by Aagaard et al. (1987) and Rudels (1987). In the still deeper layers, about 0.54 Sv of saline Arctic water flows out of Fram Strait, which is a little lower than that estimated by Heinze et al. (1990).
Note that without climate restoring a prognostic model has to adjust to all the conditions that affect the ice–ocean system, such as the salt and heat fluxes at the surface and the lateral interfaces. Any excess or deficit in fluxes of heat or salt is likely to induce a runaway system that drifts rapidly away from the climatological state. However, the standard model used here proves to be able to maintain a reasonable balance of salt and heat in the Arctic over an integration period of seven decades, and the oceanic state stays close to the climatology. This indicates that the model adapts itself reasonably well to the surface and lateral freshwater fluxes, such as precipitation, river runoff, ice growing and melting, and inflow and outflow at the ocean openings. In addition, the oceanic processes simulated by the model transport oceanic heat and regulate the heat flux at the surface and lateral boundaries in such a way that the heat balance is roughly maintained both in the upper ocean and in the deeper ocean, with only a slight heat gain in the deeper ocean.
Overall, the prognostic model behaves reasonably for many aspects of the ice–ocean circulation. It does, however, have its problems. The most notable one is that it generates unrealistically high saline water in the central Arctic that changes the density structure and affects the ocean circulation. Although the predicted mean ice growth rate is close to the estimation by Maykut (1982), ridged thicker ice, which is often in a state of melting even in winter, is not taken into account in the two-level ice model and may not be estimated by Maykut either. The absence of ridged ice, which may account for more than half of the total Arctic ice volume (Flato and Hibler 1995), may result in more ice growth than in reality, and, in turn, in more salt rejection from the ice. The salty water in the central Arctic does not seem to be properly channeled away by the oceanic circulation or to be diluted by precipitation or fresher Pacific water from Bering Strait. This is probably the reason for the high salt concentration in the central Arctic. This scenario underscores the difficulty of long-term prognostic simulation without climate corrections and the necessity of using realistic surface fluxes. To reproduce realistic surface fluxes, it is necessary to improve the system’s thermodynamics, which requires accurate representation of precipitation and radiative and turbulent surface fluxes. Model improvements are also necessary. Using a multicategory ice thickness distribution model (Flato and Hibler 1995) with ridged ice production instead of the present model with two thickness categories, for example, may improve the ice growth calculations.
Compared with the standard model, the different salinity- and temperature-restoring schemes used in the three major diagnostic models appear to have a significant impact on the predicted Arctic ice–ocean circulation. The impact comes from a chain reaction. First, the climate restoring terms change the oceanic salinity and thermal state in the mixed layer and below. The changed density structure, in turn, affects almost every aspect of the ocean circulation by altering the ocean’s dynamics and thermodynamics, which results in different surface currents and deep-water circulation and changes the oceanic processes that play a role in heat and mass transport. The ocean circulation then considerably affects the ice motion, ice outflow, and, particularly, ice thickness by altering the ice dynamics and thermodynamics. The study shows that the condition of the upper ocean and the condition of the deeper ocean, particularly the intermediate layers, are closely related despite different timescales. Altering either of these conditions will definitely affect the other and eventually affect the ice conditions in a coupled system. Consequently, Arctic modelers need to focus particular attention on the fluxes at the surface and on the oceanic processes in the intermediate layers to obtain more realistic ice–ocean circulation.
The results of the models with climate restoring particularly stress the importance of convective overturning and vertical and lateral advection to the ice–ocean thermodynamics in the intermediate layers of the Arctic and adjacent oceans. If the simulated mixed layer is excessively fresh all year long, because of strong restoring to the Levitus climatological data or for other reasons, less overturning occurs. This tends to trap more heat in the intermediate layers. Further, if the ocean velocity field is such that the vertical advection does not deliver appropriate oceanic heat upward, the ocean tends to keep even more heat in the intermediate layers. These warmer intermediate layers tend to bring more heat into the Arctic by lateral advection at Fram Strait owing to the Atlantic water inflow in these layers. With the warming up of these layers, more heat will be released from the ocean by the overturning process because of its negative feedback mechanism. However, if the heat brought in by lateral advection is not properly released by overturning and other vertical oceanic processes, the intermediate layers of the Arctic Ocean become overly warm. The heat contained in these warm layers is increasingly released and, added to the increasing lateral heat flux, gradually increases ice melting, resulting in a thinner ice cover. We speculate that, in a longer integration than the 70-yr period, the predicted temperature of the interior ocean might continue to climb, and the ice pack might continue to shrink. Thus, maintaining a proper level of convective overturning and achieving a realistic ocean circulation that predicts adequate vertical advection in the North Atlantic and in the Fram Strait area are crucial to maintaining realistic Arctic ice–ocean conditions and preventing a rapid climate drift. Otherwise too much heat can enter the Arctic, producing too warm an Arctic interior and too little ice extent, which may be an important implication for global climate models that include the prediction of the polar ice–ocean systems.
Contrary to intuition, introducing salinity restoring into the models does not seem to bring the ocean’s thermal state close to the climatology; instead it interferes with the ocean’s self-adjustment to surface and lateral flux conditions. One sensitivity study (ML-S) creates an ocean with rapidly increasing temperature, and the other (Dp-S) creates one with decreasing temperature. It is expected that if a much longer integration is performed, these two models would drift still further away from the initial climatological state. To avoid that, it is necessary to also restore the temperature below the mixed layer (as done by Dp-ST) so that the oceanic thermal state will rapidly approach a steady state and stay close to the climatology. In addition, the ocean circulation predicted by Dp-ST seems to be more coupled to bathymetry and looks more like the “observed” or climatological circulation. Therefore, for a numerical study subject to considerable limitation in computer power and interested only in short-term variations from the climatological mean field, weakly restoring (with a 5-yr restoring constant) both the ocean’s salinity and temperature below the mixed layer to climatological data is desirable because it helps the model reach a steady state quickly in both temperature and salinity while revealing short-term variations in the upper ocean. However, be aware that the oceanic processes of convective overturning and vertical diffusion may have peculiar behaviors, as shown in Fig. 12. Constraining both temperature and salinity to climatology in the mixed layer (ML-ST) does not help to bring the ocean heat budget into balance in a relatively short time because it is the temperature in the deeper layers that tends to drift away. Although the results are not shown in this paper, the ice–ocean circulation and oceanic heat transport predicted by ML-ST behave in the same way as those predicted by ML-S, only the interior ocean temperature climbs even higher and the ice pack becomes smaller.
Funding for this work was provided by NSF Grant OPP-9203470 and NASA Grant NAGW 2407. We thank J. E. Walsh and W. L. Chapman for kindly providing the sea-level-pressure and surface-air-temperature data; P. Ranelli for his help on the numerical modeling procedure; C. F. Ip for providing a program for buoy data analysis; J. Waugh, A. Schweiger, and I. Rigor for computer assistance; and G. Holloway for constructive comments on the manuscript. Finally, we are grateful to the anonymous reviewers, whose comments considerably improved this study. The plots were made using IDL and NCAR graphics.
Embedding of Mixed Layer Model
The mixed layer model is based on that of Kraus and Turner (1967), which was systematically presented by Niiler and Kraus (1977); additional information can be found in Houssais (1988) and Lemke and Manley (1984). The integrated one-dimensional mixed layer conservation equations for temperature and salinity are written as follows:
where Tm and Sm are the mixed layer temperature and salinity, D is the mixed layer depth, and −D is the location immediately below the mixed layer base. The minimum mixed layer depth is set to be 10 m, which equals the thickness of the first ocean level (see Table 1). The first terms in the right-hand sides of (A1) and (A2) are the heat and salt fluxes, respectively, at the surface and the second terms are the fluxes at the base of the mixed layer.
The surface heat and salt fluxes are determined from the forcing at the atmosphere–ocean or ice–ocean interface. The heat flux is calculated in the same way as by Hibler and Bryan (1987), whereas the calculation of the salt flux is slightly different and is written as follows:
where ρI is ice density, P is precipitation, E is evaporation at the surface, and SI is ice salinity. Ice salinity depends on ice thickness, ice type, and time of year. Cox and Weeks (1974) estimate that SI ranges from 0 to about 14 ppt (parts per thousand) for cold, level, first-year ice. Since the ice model used here does not distinguish between ice types, SI is chosen to be 7 ppt for simplicity. Note that the freshwater input due to river runoff is not expressed in (A3), it is treated as a local negative salt flux.
The fluxes at the bottom of the mixed layer are expressed as
where we is entrainment velocity determined by
where b = −g(ρ − ρr)/ρr = g[α(T − Tr) − β(S − Sr)] is the buoyancy, with ρr, Tr, and Sr being reference density, temperature, and salinity, respectively; α is the thermal expansion coefficient; β is the expansion coefficient for salinity; Δb is the buoyancy jump at the mixed layer base; u* = (|τa + F|/ρ)1/2 is the friction velocity; (b′w′)0 = g[α(w′T′)0 − β(w′S′)0] is the surface buoyancy flux; m and n are coefficients of decay of the energy production; and q2 is the rate of energy needed to agitate the entrained water. According to Kim (1976), q2 is determined by q2 = 9 max(10−4 m2 s2, u*2). According to Lemke (1987), there is a correction to the buoyancy jump for an ice-covered ocean in taking into account the effect of ocean heat flux, so the complete buoyancy jump is expressed as
where LI is the latent heat of melting of ice and c is the specific heat of seawater.
Note that the energy decay coefficients n and m in (A5) are important in that they control the rate of energy or buoyancy dissipation in the upper ocean and, hence, the extent of vertical mixing. Following the method of Lemke and Manley (1984), the dissipation process is parameterized by the following equations to ensure an exponential decay of the energy production as mixed layer depth D increases: n = exp(−D/dc), m = exp (−D/dw). A wide range of values for constants dc and dw have been seen in a number of studies (Lemke and Manley 1984; Lemke 1987; Houssais 1988; Houssais and Hibler 1993). Particularly for the Arctic Basin, Lemke and Manley tested a number of values, ranging from about 9 m to less than 30 m; for the Greenland Sea, Houssais used larger values. For simplicity, dc and dw are uniquely expressed here as linear functions of D regardless of region:
These relationships are tuned to obtain a reasonable mixed layer in the Arctic that does not become too deep in winter.
For a one-dimensional mixed layer model, (A5) can be used to predict the mixed layer depth during entrainment (we > 0) by dD/dt = we. During detrainment (we ⩽ 0), a new equilibrium mixed layer depth is determined diagnostically by solving (A5) with we = 0. For a three-dimensional mixed layer model, the effects of horizontal advection have to be taken into consideration to determine the actual mixed layer depth:
where um and υm are the horizontal velocity components in the mixed layer.
The basic objective of the embedding procedure is to couple the vertical mixing calculated by the mixed layer model with the advective and diffusive processes simulated in the multilevel primitive equation model. Since the depth of the mixed layer does not necessarily coincide with any of the fixed ocean levels, work has to be done to ensure that incorporating the steplike vertical profile of the mixed layer density will not cause loss of the conservation of heat, salt, and momentum, both in the mixed layer model and in the ocean model, during entrainment or detrainment. The steps in the coupling procedure are outlined below:
1) Run the ocean model to obtain U, T, and S due to advective and diffusive processes at every level of the ocean model.
2) Compute the contribution of the advective changes to the mixed layer depth using (A7) without entrainment we.
3) Adjust the mixed layer temperatures and salinities obtained in step 1 to the mixed layer depth obtained at step 2 (Resnyanskiy and Zelen’ko 1991):
where k is the level containing the mixed layer base, as shown in Fig. A1, and Tl and Sl are temperature and salinity at level l in the ocean model.
4) Correct the ocean model’s profiles of T and S to ensure the conservation of the related quantities:
5) Calculate the mixing process: If we > 0 from the solution of (A5), new mixed layer is defined using (A7) without the advection and diffusion terms, which have already been used at step 2. If we ⩽ 0, the mixed layer normally retreats to a shallower depth. In this case, the embedding technique of Adamec et al. (1981) is used, which adjusts Tm and Sm based on additional potential energy conservation (see Adamec et al. 1981 for details). This adjustment may work well in reforming a new, shallower mixed layer in equatorial or temperate oceans, where the temperature of the upper ocean is usually well above the freezing point. However, when the new mixed layer becomes much shallower than before, say, compared with a previous time step, a problem arises with the adjustment for ice-covered oceans or areas that have a mixed layer temperature close to freezing; the adjustment may result in such a temperature decrease that the mixed layer temperature falls below freezing, causing a phase change and possible formation of a large volume of ice that is not realistic. To avoid this problem, Houssais and Hibler (1993) dropped the requirement of potential energy conservation, which seems to work well in the subarctic ocean. Here Adamec et al.’s technique, with potential energy conservation, is still used, but every time the diagnostic mixed layer depth is obtained by solving (A5), the average of that depth and the one obtained in the previous time step is used as the new depth. This was found to greatly limit the possible dramatic temperature decrease in the mixed layer while delaying the mixed layer approaching its shallowest position by only a few time steps.
Imposition of Open Ocean Boundary Conditions
The following measures are taken to introduce open boundary conditions for the Bering and Denmark straits and Faroe–Shetland Passage:
1) 30-day relaxation is imposed at open boundaries and in their vicinity.
2) ∂T/∂n = ∂S/∂n = 0, where n is the direction perpendicular to open boundaries.
3) The vertically integrated transport streamfunctions (ψ) that are specified along the open boundaries are shown in Table B1. The vertically averaged velocity, or external mode of velocity (see Semtner 1986), at the ocean boundaries can be determined from the specified streamfunction:
where H is the ocean depth.
4) A simplified version of the internal mode of velocity (i.e., the velocity deviation from the external mode) can be obtained from the following geostrophic equations:
where ps is the pressure at the ocean surface and p0 is hydrostatic pressure component. Let
at the open boundaries.
Corresponding author address: Dr. Jinlun Zhang, Applied Physics Laboratory, University of Washington, 1013 NE 40th Street, Seattle, WA 98105-6698.