[June 14, 2013] Effects of Stochastic Ice Strength Perturbation on Arctic Finite Element Sea Ice Modeling [Journal of Climate] Tweet (Journal of Climate Via Acquire Media NewsEdge) ABSTRACT The ice strength parameter P* is a key parameter in dynamic/thermodynamic sea ice models that cannot be measured directly. Stochastically perturbing P* in the Finite Element Sea Ice-Ocean Model (FESOM) of the Alfred Wegener Institute aims at investigating the effect of uncertainty pertaining to this parameterization. Three different approaches using symmetric perturbations have been applied: 1) reassignment of uncorrelated noise fields to perturb P* at every grid point, 2) aMarkov chain time correlation, and 3) aMarkov chain time correlation with some spatial correlation between nodes. Despite symmetric perturbations, results show an increase of Arctic sea ice volume and a decrease of Arctic sea ice area for all three approaches. In particular, the introduction of spatial correlation leads to a substantial increase in sea ice volume and mean thickness. The strongest response can be seen for multiyear ice north of the Greenland coast. An ensemble of eight perturbed simulations generates a spread in the multiyear ice comparable to the interannual variability of the model. Results cannot be reproduced by a simple constant global modification of P*. (ProQuest: ... denotes formulae omitted.) 1. Introduction Because of the diversity in spatial and temporal scales of physical processes and the limited resolution in climate models, parameterizations are generally needed to describe the influence of small-scale processes on largescale flows. This influence is mostly considered in terms of the mean of the models' prognostic variables, while the effect of unresolved scales on higher-order moments such as the variance is frequently ignored. Because of the numerical implementation of the parameterizations, model errors arise. Especially in numerical weather prediction (NWP) the idea of applying probabilisticmethods to previously deterministic models to address such model errors and to increase the model spread has become a well-established research topic in recent years. As summarized by Palmer (2012), the inclusion of stochastic parameterizations may help to better estimate model uncertainty and reduce biases or even the uncertainties themselves. In this context, Buizza et al. (1999) applied a stochastic perturbation to the parameterized tendencies in their European Centre for Medium-Range Weather Forecasts (ECMWF) Ensemble Prediction System(EPS) to address model error within the parameterizations of physical processes. This successfully increased the ensemble spread. This was followed by amore sophisticated spectral stochastic kinetic energy backscatter scheme (SSBS) (Berner et al. 2009), based on stochastic perturbations of streamfunctions. Another approach that also aims at the inclusion of model error within the physical processes in NWP has been applied by Lin and Neelin (2000), Bright and Mullen (2002), and Li et al. (2008). These authors stochastically perturbed parameters of a certain parameterization scheme to account for uncertainties that are related to a deterministic parameter choice, as some of the parameters do not have an obvious analog in nature. A first-order Markov process is used by Lin and Neelin (2000) to perturb the convective available potential energy within the parameterization of deep convection. The same concept of temporal correlation is applied by Bright and Mullen (2002) to perturb the grid-scale vertical velocity at the lifting condensation level within the cumulus parameterization and the critical bulk Richardson number within the planetary boundary layer scheme. Li et al. (2008) employed a more sophisticated parameter perturbation in their convective and condensation schemes. They applied a first-order Markov process to create randomnumbers as spectral coefficients of spherical and Fourier harmonics expansions in the horizontal and vertical direction, respectively. Thereby, spatial correlation is introduced to the stochastic parameter perturbations. Sea ice models contain several sensitive parameters that need to be chosen carefully and are used to tune the model. Early studies-for example, by Stössel et al. (1990), Owens and Lemke (1990), and Harder (1996)- show that the ice strength parameter P* of the viscous- plastic sea ice rheology introduced by Hibler (1979) is among the important parameters within the sea ice dynamics, together with parameters such as the drag coefficients for the momentum fluxes between the ice and the ocean and atmosphere. Sensitivity experiments carried out prior to this study have confirmed the importance of P* when it comes to sea ice distribution and driftpatterns in our model. The quantity P* determines the internal ice strength; for higher values of P* the ice is less easily deformed in case of a convergent drift. Changes in the sea ice driftinfluence the advection of sea ice, which has an impact on the ice thickness distribution, especially on the formation of pressure ridges that often accrete to multiyear ice, and on the ice concentration. As sea ice plays an important role for the surface albedo and has a strong effect on the transport of heat and momentum between atmosphere and ocean, changes in the ice coverage influence the thermodynamics of the model and in coupled models affect both the ocean and the atmosphere components. So, the parameter P* has an impact not only on the sea ice driftitself and on the thickness distribution, but also on aspects such as the Arctic freshwater export linked to sea ice driftthrough the straits, the heat budget of the ocean, sea ice production rates, and bottom water formation. As a consequence, a proper representation of the sea ice dynamics is necessary for accurate predictions of climate in the high latitudes. The value of P* is only weakly constrained, however, as it cannot be measured directly and may vary strongly in time and space (Harder and Fischer 1999). Tremblay and Hakakian (2006) estimated upper and lower bounds for P* based upon the deformation law of Hibler (1979) using satellite data for sea ice driftand reanalysis data for sea level pressure. Even though they found values of the same magnitude as other studies that used in situ measurements or driftbuoy data, the derived values vary quite considerably. They also state that the timeaveraged atmospheric forcing used in any given model as well as the choice of parameters such as the air-ice drag coefficient may have a strong influence on the appropriate value of P*. This leads to further uncertainties. Because of these uncertainties and the global averaging of a highly variable parameter, P* seems a good first candidate for applying a stochastic approach as a method to tackle the model error connected to the physical process of subgrid-scale sea ice deformation and to investigate the changes arising in the sea ice thickness distribution due to the nonlinearity of the parameterization. In this work we present a first approach to include uncertainty estimates in a sea ice model by the use of stochastic parameterizations. We introduce a new suite of schemes to stochastically perturb the ice strength parameter P* and demonstrate that this symmetric parameter perturbation leads to an increase in total Arctic sea ice volume and a distinct change in the local Arctic sea ice thickness distribution. Adding first temporal and then spatial correlation to the perturbations leads to an increasing impact on the Arctic sea ice coverage. The spatial correlation in particular will prove to be of great importance in the perturbation scheme. While the influence of the seasonal cycle reduces the effects of the perturbed sea ice strength parameterization for firstyear ice and along the ice edge, the multiyear ice exhibits a strong sensitivity to the new scheme and tends to increase considerably in thickness. The Arctic monthly mean sea ice thickness in an ensemble simulation shows an increase of about 15% in January and about 111% in September compared to the reference simulation. The spread in multiyear ice thickness generated by the parameter perturbation is comparable to the interannual variability of the model and may give an indication for a range of values for P* to be used in data assimilation studies. The range chosen for most perturbations in this study is P* 2 (5000Nm22, 35000Nm22). This paper has the following structure. In section 2 we will give a short overview of the Finite Element Sea Ice- OceanModel (FESOM) (Wang et al. 2008; Timmermann et al. 2009; Sidorenko et al. 2011) of the Alfred Wegener Institute for Polar and Marine Research (AWI) and the sea ice rheology that is used by the model. Section 3 describes the three different parameter perturbation approaches that have been studied: one with no spatial correlation between any two nodes and just a simple form of temporal correlation, a second one again with no spatial correlation but a temporal correlation given by a firstorderMarkov process, and finally a third version featuring the temporal correlation according to the Markov process of the second approach and also a spatial correlation between nodes with common neighbor nodes. Results of the three approaches in comparison with a deterministic reference simulation and with observations of sea ice draftand concentration will be analyzed in section 4, focusing on the Arctic Ocean. Finally, section 5 will provide conclusions and a short outlook for future work. 2. Model background We use the Finite Element Sea Ice-Ocean Model to perform simulations of the Arctic sea ice. Spatial discretization of the model is implemented via finite elements on a triangular grid at the ocean surface and a tetrahedral grid with z levels for the three-dimensional ocean modeling. The ocean model discretizes the momentum, heat, salt, and mass continuity equations based on the hydrostatic assumption and is coupled to a sea ice model. The model setup in this study is forced by the interannual varying Common Ocean Reference Experiment (CORE), version 2 atmospheric dataset by Large and Yeager (2009). FESOM under the normalized year forcing of CORE I has been evaluated and compared to other models by Sidorenko et al. (2011). Further details on FESOM have been presented by Wang et al. (2008) and Timmermann et al. (2009). In the sea ice model the prognostic variables effective (mean) sea ice thickness hice, lateral sea ice velocities uice and yice, and sea ice concentration A, as well as mean snow layer thickness hs, are calculated. For this purpose, prognostic equations for hice, hs, and A describing the rate of change due to freezing and melting processes and advection are solved. The dynamical part of the model describes the rate of change in the sea ice velocities via the momentum balance (Timmermann et al. 2009). Sea ice and ocean models use the same surface mesh, which makes their coupling straightforward. The global finite element mesh is locally refined in the Arctic and especially along the coastlines and in straits, the resolution ranging from over 300km in the Pacific Ocean to about 11km at the coasts. The mesh of the Arctic region can be seen in Fig. 1. The time step in this study is Dt 5 1h. The momentum equation for the sea ice velocities discretized in FESOM is ... (1) withmbeing the mass per unit area, ui the horizontal sea ice velocity vector, tair and tocean the atmospheric and oceanic stress, respectively, f the Coriolis parameter, k5(0, 0, 1), g the gravitational acceleration, and ho the sea surface height of the ocean. For this study the important term in Eq. (1) is the last term on the right-hand side that describes the internal forces Fint 5(F1, F2) within the sea ice. For the computation of the internal stress FESOM offers the classical viscous-plastic rheology introduced by Hibler (1979) with adjustments by Harder (1996), and alternatively the elastic-viscous-plastic rheology devised by Hunke and Dukowicz (1997) following in its implementation the documentation of the Community Ice Code (CICE) model (also known as the Los Alamos Sea Ice Model; Hunke and Lipscomb 2010). The latter has been used for this study. Both rheologies describe the internal forces as the divergence of the two-dimensional stress tensor s: ... (2) For the elastic-viscous-plastic rheology the tensor components sij are described by the differential equations ... (3) where e is the deformation rate tensor ... (4) h and z are the nonlinear shear and bulk viscosities, respectively, dij is the Kronecker symbol, and E is some elasticity constant (Hunke and Dukowicz 1997). The ice strength P is given by ... (5) ... (6) where Dmin represents a small regularization value that creates a smooth transition from viscous to plastic regimes, following the idea of Harder (1996), and e is the eccentricity of the elliptical flow curve (i.e., it describes the relation between z and h 5 z/e2). The important parameterization concerning the perturbations introduced in the next section is that for the sea ice strength given by ... (7) with the dynamic ice strength parameters P* and C. The viscosities h and z also depend linearly on Pp. Equations (5) and (7), through (2) and (3), describe the resistance of the sea ice against deformation. For small deformation rates the internal forces increase with rising deformation rates, in a manner that depends on the value given by Eq. (7). But as soon as a limiting value of D is reached, a further increase in deformation rates does not increase the internal forces. The sea ice behaves plastically and, by way of its drift, tends to pile up. Therefore Eq. (7) influences the behavior of sea ice especially in convergent drifts, where hice and A have an impact on the rate at which the sea ice is deformed and thereby increases in thickness. This highly nonlinear dependence of internal stresses on ice strength and hence P* will be analyzed in the following section by a symmetric parameter perturbation. 3. Parameter perturbation The stochastic perturbation in this work is applied to the sea ice strength parameter P* governing the internal forces within the sea ice dynamics. In the purely deterministic version of FESOM the parameter P* is fixed in time and space at a value of P*5Pref * 520 000Nm22. This does not account for small-scale variations imposed by, for example, different local melting and freezing histories or brine pocket formation. Furthermore, the associated natural variability in space and time around a givenmean in sea ice thickness and compactness is not resolved. Some of this variability can be introduced by a symmetric perturbation of P* around Pref * . The general approach is based on perturbing Pref * via ... (8) The indices i and j give the nodal index in space and the time step, respectively. The value x(i, j) is a random value from a preassigned distribution, bounded within some range 2a , x(i, j) , a, with 0 , a # 1. For the use in the sea ice rheology in FESOM, nodal values are averaged over elements where rheology calculations take place. This general kind of parameter perturbation has been used in different parameterizations within atmospheric general circulation models, for example by Lin and Neelin (2000) in a convective parameterization or by Bright and Mullen (2002) in cumulus and in planetary boundary layer parameterizations. In both studies the random numbers x(i, j) have been correlated in time, but not in space. In this work we will present three different possibilities for the parameter perturbation in Eq. (8). In the first approach, reassignment of random values will occur after a predefined time interval and there will be no correlation between different nodes. For the second approach the reassignment of random values will be replaced by a time correlation via a Markov process. Finally, the third approach will include a limited correlation in space as well as the correlation in time via a Markov process. In this context, the inclusion of correlation in time and space seems reasonable. Sea ice has a strong memory when it comes to processes such as formation, melting, and deformation, all of which have an impact on the resistance of the ice toward further deformation.And even though small-scale structures such as ridges and leads occur in the ice, there is generally considerable correlation within the spatial ice thickness distribution, which is connected to different ice strength regimes. The relevant parameters of spatial and temporal correlation for the new stochastic schemes will receive their values from references given in section 4. a. Reassignment time step perturbation (RP) The RP approach is implemented as follows. At the end of each time interval with length (Dt)pert, every ice-covered node within the model area is reassigned a random number from some distribution. In addition, open ocean nodes that freeze over between those larger time steps receive some random number as well. Gaussian distributions are widely used for parameter perturbations (Lin and Neelin 2000; Bright and Mullen 2002). The interval for possible values of the random numbers is bounded, however, because P* has to stay within physically realistic limits. Therefore the traditional Gaussian distribution cannot be used. In this study three differentGaussian-like distributions have been employed. Appendix A gives a short overview over the so-called x- and y-truncated Gaussian distributions (Bardsley 2006) and a third distribution, which is generated through exponential transformation of a Gaussian distribution. The important parameters for this first approach are the limit a of the domain of the random number, the reassignment time step (Dt)pert, the truncated or bounded distribution function used, and the variance of the underlying Gaussian distribution function. Values will be assigned in section 4. b. Markov process time correlated perturbation (MTP) To include some correlation in time that does not only depend on a reassignment time step, the approach of section 3a is adjusted. Following Bright and Mullen (2002) and Lin and Neelin (2000) a Markov process of the general form ... (9) is used to keep correlation of the random value y(i, j) at the current time step with the previous one, y(i, j 2 1), with y(i, 0)50 as initial condition and also y(i1, j121)50 if j1 is the first ice-covered time step for node i1 after an ice-free period. In Eq. (9) 0 , a , 1 denotes thememory of the first-order Markov process (i.e., the autocorrelation). As derived in appendix B, ... (10) where l 5 1/t and t is the relaxation time. The perturbation z(i, j) is an independent Gaussian random number with zero mean. We set the variance of the z(i, j) to ... (11) with s being comparable to the standard deviation of the Gaussian perturbation used in section 3a, (Dt)pert the reassignment time step also used in section 3a, and Dt themodel time step. Contrary to the approach in section 3a, the perturbations z(i, j) are applied at each time step. To use the random values produced by Eq. (9) the values need to be transformed into the respective values x(i, j) with limited range. How this is done is also described in appendix B, for all the three possible bounded distributions. Differences in the results of the runs due to the three different corresponding transformations will be discussed in section 5. c. Markov process time and space correlated perturbation (MTSP) To include some spatial correlation in the perturbation of P*, the random number z(i, j) inEq. (9) is adjusted as ... (12) Here j(i, j) and j(k, j) are N(0, 1) distributed random numbers, q(...) is a weighting function depending on the distance dik in meters between the current node i and its neighbor k, and Ni is the set of indices of the direct neighbor nodes of node i. Dividing by the square root of the sum in (12) normalizes the variance of the random number z(i, j) so that the corresponding random variable is still Gaussian distributed with the variance (11) and zero mean. The weights (...) are exponential functions given by ... (13) with some correlation distance dcorr in meters. In this approach the variance of a random variable that corresponds to (12) is given by Eq. (11), while the covariance between two nodes is greater zero if they have common neighbors. Otherwise the covariance is equal to zero. The procedure for transforming the values y(i, j) into the respective values x(i, j) with limited range is the same as for MTP as all variance calculations given in appendix B stay the same. 4. Model setup The three different approaches of section 3 for the parameter perturbations have been implemented in the framework of FESOM. To analyze the changes caused by the perturbations, a reference simulation REF with fixed parameter P*5Pref * 520 000Nm22 has been spun up for 12 yr. It was continued further and also used to initialize the different test cases for the stochastic perturbations. For the first part of this study we compared the total Arctic sea ice volume and area as well as the monthly mean sea ice thickness and concentration distribution of one realization of each stochastically perturbed model setup with the reference run. Integration time for these sensitivity experiments was 3 yr, starting from 1990. As will become clear in section 5a, the approach with the strongest impact is MTSP, with Markov process time and space correlated perturbations, which is the reason why the sensitivity studies presented here will focus on this last approach. For the other two approaches, RP and MTP, only one experiment each will be discussed in section 5a. Table 1 shows the different model setups, introduces their names, and summarizes their impact on total Arctic sea ice volume and area. Parameter ranges for the relaxation time t follow Lemke et al. (1980), where the general idea of the Markov chain has been used to estimate, among other things, relaxation times for sea ice extent anomalies in the Arctic and Antarctic. The lower values of 2 h and 1 day have been tested as well to reduce the temporal correlation considerably. The range for correlation distances follows estimates of spatial scales for a valid interpretation of the central Arctic sea ice thickness distribution function suggested by Flato (1998). The parameter range ((12a)Pref * , (11a)Pref * ) for possible values of P* follows several publications concerning sensitivity studies of P*, among other parameters, within different models (Owens and Lemke 1990; Harder and Fischer 1999). For this study we used from P* 2 (10 000Nm22, 30 000Nm22) to P* 2 (3000Nm22, 37 000Nm22). The value for the variance s2 is matched to the value of a to avoid the danger that the distribution function of the random values becomes too flat after transformation and therefore tends to lose its Gaussianlike shape. Values of s that have been tested but will not be discussed here are 0.25 and 0.5, as similar effects can be generated by changes in other parameters [e.g., (Dt)pert]. The last parameter that needs to be chosen is (Dt)pert. We chose values for (Dt)pert to be within the range of 1 to 7 days to simulate the differences between the accumulated influence of a diurnal up to a weekly perturbation cycle. Also a value of 1 h has been tested. For the analysis of an ensemble run in section 5b we used an ensemble ENS of eight members with the configurations given in Table 1 for MTSP0. Integration time for the ensemble was 17 yr, also starting from 1990. In addition, we ran three deterministic experiments with fixed parameter values for P* (REF1 with 12 500Nm22, REF2 with 15 000Nm22, and REF3 with 17 500Nm22) for 6 yr, to be able to compare the changes caused by the parameter perturbations in the ensemble with the effect of a uniform reduction of P*. 5. Results a. Sensitivity studies In Fig. 2 results from experiments RP0, MTP0, and MTSP0 are shown, compared to the change caused by a global reassignment of P* 5 15 000Nm22 in REF2. The differences in monthly mean sea ice volume and concentration between each approach and the reference run REF for the Northern Hemisphere for 3 yr reveal that the time and space correlated perturbations in MTSP0 are by far the most influential, comparable in magnitude to the deterministic REF2 experiment with P* 5 15 000Nm22. It can also be seen that all three approaches lead to an increase in total ice volume and a decrease in total ice area. This is consistent with the reasoning by Owens and Lemke (1990): If P* values in the stochastic simulations are low, sea ice thickness increases locally along lines of convergent driftwhile the sea ice coverage around those areas is decreased. This allows ice-free areas to freeze over and the cycle is repeated. Locally high P* values, on the other hand, do not enforce a decrease in sea ice thickness within that area but simply reduce convergent sea ice driftand keep the ice thickness at the current level. Therefore, even for symmetric perturbations, the tendency is toward locally increased sea ice thickness and decreased sea ice concentration, which leads to a slight increase in total Arctic sea ice volume. For the cases of global parameter changes of P* (e.g., REF2) all Arctic sea ice is weakened compared to the reference run REF. This leads to the growth of thicker ice through convergence and a resulting opening of sea ice cover. While this is a large-scale effect, in the simulations with symmetric perturbations of P* the effects are very localized, which might lead to a hindered increase in sea ice thickness at one point in time and space (high P* value) but also to a more than compensating increase in ice thickness at another point (low P* value). Generally the effect of the stochastic perturbations on the sea ice volume is largest during the melting season in spring, while the total sea ice area shows the strongest response in the summer months, as will be discussed in section 5b for the ensemble simulation. For MTSP0 Fig. 3 shows the evolution of the mean monthly difference in total ice (right) area and (left) volume with respect to the reference run REF for 17 yr. It takes about 3 yr until the perturbed simulation run reaches a new quasi-steady state. There is an anomaly in the modeled sea ice volume difference around the year 1997 due to a rapid decrease in the total REF and also MTSP0 sea ice volume in 1995 and 1996. The mean total change in Arctic ice volume for the 17 yr simulated (10.83 3 103km3) makes up about 4% of the reference run's sea ice volume. For the total sea ice area the mean value of the decrease is 20.12 3 106km2, which is only about 1% of the reference run's mean sea ice area. Figure 4 shows the local effect of the three different approaches RP0, MTP0, and MTSP0, revealed by the differences in monthly mean sea ice thickness and concentration between the stochastic model runs and the reference run for the summer 1992 [June-August (JJA)]. As can be clearly seen, the magnitude of the local thickness and concentration changes is increased by the inclusion of time and especially space correlation within the stochastic perturbations. Especially along coastlines, in areas of convergent driftand withmultiyear ice coverage, the sea ice thickness is increased. As the sea ice concentration is high in the centralArctic Ocean, changes of ice concentration mostly occur along the ice edge and in areas of lower concentration where there is also convergent drift(increase of concentration) or even divergent driftor reduced convergent driftrelative to the reference run, newly created by driftpattern changes in nearby areas. For winter months the changes in sea ice concentration are mostly confined to a small front parallel to the eastern Greenland coast and within the Bering Strait, as modeled sea ice concentrations farther north are around 95%-100% in winter (not shown) and thus do not leavemuch space for a further increase of concentration in convergent drift. In summary, the effect of the stochastic perturbations on the local sea ice concentration is rather small, with the strongest influence during the summer months. The influence on the location of the ice edge is very limited, as the changes in sea ice concentration especially during the winter months but also in summer are for most places not strong enough to fully remove the sea ice cover in the monthly means. The choice of the distribution function used for stochastic perturbations only has a moderate effect on changes in total ice concentration and volume. Especially for MTP and MTSP, the transformation into the x-truncated Gaussian distribution shows the biggest difference in comparison to the reference run (not shown). This is because the truncated Gaussian distribution function assumes a relatively flat shape with the time-step-dependent increase in variance as given by Eq. (B3) (see appendix B). As mentioned before, the small P* values have a bigger influence on the changes in sea ice concentration and volume than the high values, both of which occur more often with random numbers that follow a distribution function with a broad maximum. From here on the x-truncated Gaussian distribution function will be generally used. Similar results with a slightly smaller difference between reference run and stochastic run can be obtained with the simpler transformation given by Eq. (A1) (see appendix A) or the transformation into the y-truncated Gaussian distribution (not shown). The last two columns of Table 1 show the mean deviations in sea ice volume and area of the stochastic experiments from the reference run REF for the time period 1990 to 1992. These sensitivity studies revealed that the most sensitive parameters in the MTSP are the limiting value a for the range of P* (MTSP1 and MTSP2) and the correlation distance dcorr especially between values of 10 to 100km (MTSP8 and MTSP0). A further increase in correlation distance to 1000km (MTSP9) does not lead to a more pronounced impact on changes in ice volume and area compared to MTSP0 with 100-km correlation distance. These results provide insight into the range of spatial correlations within the perturbations that are influencing the general sea ice thickness and concentration patterns. The correlation time t for MTP and MTSP has a quite small effect when it is varied in the limits investigated by Lemke et al. (1980). Only very small values can decrease the influence of the temporal correlation in a way that reduces the effect toward a no-time-correlation approach (e.g., MTSP5 and MTSP6). Changes in the redistribution time interval (Dt)pert forRPalsohave a small influenceon the difference between reference run and stochastic run, with an increasing difference occurring with an increasing redistribution time interval size (not shown). The parameter that has a strong effect on the form of the distribution function after transformation is the variance ~s2 of the applied perturbations in MTP and MTSP, which itself can be regulated through the parameters (Dt)pert (MTSP3 and MTSP4) and s [see Eq. (11)]. As the time dependent total variance given by (B3) rises quickly in the beginning but then tends asymptotically toward its limiting value (B4), the limiting value is of importance when it comes to the form of the distribution function that most long-time ice-covered nodes share after transformation. Most values for ~s used in this study are relatively high. Because of the design of the transformations, the resultant form of the chosen bounded distribution function does not vary much for high values of the total variance (B4) of the underlying Gaussian distribution. In our sensitivity studies medium to high values for variance lead to the biggest difference in sea ice volume and extent. Very high total variance values that are for example reached with (Dt)pert 5 1h (MTSP3) do not only seem unlikely but also do not show a stronger impact than values that are quite a bit smaller because the transformed distribution functions do not differmuch. Only small ~s and therefore total variance values actually lead to a considerable decrease in the impact of the stochastic perturbations [e.g., with (Dt)pert 5 7 days (MTSP4)]. b. Ensemble run 1) GENERAL EVALUATION The ensemble experiment with eight ensemble members is used to analyze the mean effect of the Markov process time and space correlated perturbations of P*. Furthermore, we are interested in the spread generated by the ensemble. The parameter choice for the members is the same as for MTSP0 in Table 1. Figure 5 shows the impact of the stochastic perturbations on the monthly mean Northern Hemisphere ensemble mean sea ice volume and area compared to the reference run for the first 6 yr. Also shown are the differences of deterministic runs with fixed parameter values for P* (REF1 with 12 500Nm22, REF2 with 15 000Nm22, and REF3 with 17 500Nm22) to the reference run REF with P* 5 20 000Nm22. It can be seen that the changes in Arctic sea ice volume (increase) and area (decrease) that occur due to the stochastic perturbations are persistent even after performing the ensemble mean and cannot be reproduced by simple recalibration of fixed parameter values. Figure 6 shows the changes in local sea ice thickness due to the stochastic perturbations for the mean of ENS (top left) compared to changes produced by the fixed parameter value run REF2 with P* 5 15 000Nm22 (top right) for March 1995. While the changes due to the constant modification of the parameter value are very smooth in their spatial distribution in the central Arctic, the stochastic perturbations act locally and thus introduce considerable smaller-scale variability. The tendency toward increased multiyear ice thickness is consistent in both experiments. This is also visible when looking at the modeled ice thickness distribution of the reference run REF (Fig. 6, bottom left), the ensemble mean (Fig. 6, bottom middle), and REF2 (Fig. 6, bottom right) for the entire Arctic.As can be seen, the ice thickness distribution of the ensemble mean shows a longer right tail similar to REF2 with P* 5 15 000Nm22 when compared with REF. Sea ice thickness in the ensemble can reach more than 9m.On the other hand, while the general structure of the distribution is the same for the three simulations, with a peak for thin first-year ice and another one for the thicker multiyear ice with a thickness of about 2m, both peaks are more pronounced for the ensemble mean (Fig. 6, bottom middle). All three distributions are reasonable when compared to measurements, asmultiyear and first-year ice are clearly distinguishable, even though the location and height of the peaks might be slightly different in measured samples depending on sampling time and place (S. Hendricks 2012, personal communication; Yu et al. 2004). The difference in monthly mean sea ice thickness averaged for the entire Arctic (monthly mean total sea ice volume divided by monthly mean sea ice area) between the mean of the ensemble run ENS and the reference run REF with Pref * ranges from about 10.08m (about 15% of reference run mean sea ice thickness) in January to more than 0.26m (more than 111%) in September, here discussed for 2004. This implies that the biggest effect of the perturbations pertains to thick multiyear ice that survives the spring melt. The ratio of the area of multiyear ice to the area of first-year ice increases in spring toward summer when the total sea ice area decreases. Moreover, the changes in multiyear ice thickness are strong enough to increase the relative difference between ensemble and deterministic run even though the mean monthly ice thickness of both experiments increases substantially during spring from about 1.7m in November to about 2.6m in July 2004. In Fig. 7 the relative influence of the ensemble run on sea ice volume and area compared to the reference run can be seen. Once again this figure shows the increased influence of stochastic perturbations in the late spring and summer months. As the effects of the stochastic perturbations in the dynamics scheme of the model add up over time when it comes to sea ice thickness, regions with all-year ice coverage produce a pronounced difference in ice thickness compared to the reference run REF. Figure 8 illustrates this for the 14-yr mean seasonal difference between the reference run REF and the mean of the ensemble ENS. The first three years (1990-92) have been omitted in the averaging for Fig. 8, as the quasi-steady state of the difference is reached around 1993. The effect of the stochastic perturbations on the multiyear ice is of comparable magnitude in winter and summer, whereas first-year ice exhibits only a minor response to the parameter perturbation in the long-term mean. In addition to the long-term means of Fig. 8, Fig. 9 presents the monthly mean sea ice thickness standard deviations of the stochastic ensemble for the Arctic for March 1996 (top left) and March 2004 (top right). Furthermore the long-term mean standard deviation within the ensemble for the months of March 1993 to 2006 is illustrated by the bottom leftpanel of Fig. 9. Its magnitude is comparable to the modeled interannual variability of the reference run, shown in the bottom right panel as the sea ice thickness standard deviation for March 1993 to 2006. The ensemble-generated spread is large in regions where dynamics has a strong influence during the whole year, while it is low in regions with a strong impact of seasonal melting and freezing processes and especially where the sea ice completely melts in summer. In the high-deformation region north of Greenland the ensemblegenerated variability exceeds the interannual variability in the reference run by a factor of 2. The small ensemble spread for regions with predominantly firstyear ice is even more visible in the summer months, even though the modeled interannual variability of REF is high in those months, especially in September (not shown). 2) COMPARISON TO OBSERVATIONS The effect of the stochastic perturbations on changes in sea ice concentration is small compared to changes in the sea ice thickness distribution. Comparing model sea ice concentration data with satellite data from Cavalieri et al. (1996), the reference FESOM model run REF tends to overestimate the sea ice concentration, especially in Fram Strait and in the Greenland Sea (not shown). The stochastic ensemble exhibits only very marginal improvements in those regions during the winter months and no noteworthy improvements during the summer months (not shown). Changes in sea ice velocities are generally small as well (not shown). The perturbations influence the sea ice velocities on rather short time scales, which has an influence on the sea ice thickness and concentration distribution, as changes concerning those variables can add up over time. Only strong changes in the sea ice thickness may (through sea ice rheology) cause changes in velocities that last for longer times. Comparing reference and ensemble experiments to satellite data for sea ice thickness [e.g., 3 October to 8 November 2004, derived by Kwok et al. (2007, 2009) from the Ice, Cloud, and Land Elevation Satellite (ICESat)], the ensemble is in better agreement with the observations for the multiyear ice north of Greenland, while the underestimation of the sea ice thickness in the eastern Arctic remains mostly unaffected (Fig. 10, bottom row). This holds for most of the five autumn and five winter months of the years 2003 to 2008, for which Kwok et al. (2007, 2009) presented observational data. The general distribution with the high sea ice thickness along the northern Greenland coast and the Canadian Arctic Archipelago is reproduced by the reference run REF as well as by the mean of the ensemble run ENS (Fig. 10, top row). When it comes to the comparison of local sea ice thickness with upward looking sonar (ULS) data, which give time series of observed sea ice draftat fixed locations, changes due to the stochastic perturbations are small and no general improvement of the ensemble mean compared to the reference run can be seen. This is mostly due to the fact that the ULS data points available for this study are located in regions where the impact of the stochastic perturbations is low (i.e., in the Beaufort Sea, Fram Strait, and Greenland Sea). As can be seen in Fig. 11 for two ULS data series, one in Fram Strait (right) and the other in the Beaufort Sea (left), the differences between the ensemble members of the stochastic model run are small and decrease in times of increased thermodynamic influence or divergent drifts, when the value of the ice strength is of little importance (e.g., around day 10 and day 130 in Fig. 11, left). The seasonal cycle has a strong effect on the ensemble spread, which goes down to zero when the sea ice melts in summer. We have performed root-mean-square and mean difference calculations between the observations [9 ULS data series in the Beaufort Sea (Melling and Riedel 2008) and 11 ULS data series in the Fram Strait and Greenland Sea (Witte and Fahrbach 2005)] and the corresponding modeled ensemble mean time series. The same has been done between observations and the modeled reference run time series. Comparing the results indicated minor improvements for some time series due to the stochastic perturbations but also some minor worsening of results for other locations. Generally the ensemble spread is very small in these regions. On the other hand, Fig. 12 shows that the sea ice thicknesses produced by the ensemble members at a location near the northern Greenland coast vary considerably. This is also true when the ensemble mean is compared to the reference run. Standard deviation of the ensemble for the time series at 848N, 408W reaches, for example, 2.25m and more in January 2002, while the mean difference between ensemble mean and reference run in, for example, the year 1994, is 0.79m (about 12% increase). The cause for this is the reduced seasonal variability produced by thermodynamic melting and freezing compared to the variability due to dynamics in the central Arctic Ocean. Furthermore, the high amount of multiyear ice is of major importance, especially for the sea ice rheology. Unfortunately, sea ice thickness data in the region north of the Greenland coast are not available in a way that allows for a meaningful comparison. This is also true when it comes to data necessary for calibrating the perturbation scheme correctly, with respect to the relation between ensemble mean error and ensemble spread. 6. Summary and outlook Three different approaches to randomly perturb the parameter P* of the ice strength parameterization of the sea ice dynamics module in FESOM have been presented, with increasing complexity concerning temporal and spatial correlation. Because of the nonlinearity of the sea ice rheology, a symmetric perturbation of P* leads to a general increase inArctic sea ice volume and a decrease in the sea ice area compared to a deterministic reference run. Inclusion of some temporal correlation via a Markov process and especially a spatial correlation between nodes with common neighbors proved to have a strong influence on the changes in the sea ice thickness and concentration distribution compared to the reference run. These changes cannot be reproduced by a global reassignment (reduction) of constant parameter values for P* in the purely deterministic model. Especially in regions where the influence of sea ice dynamics is high and where there is no complete melting of sea ice during summer, the most sophisticated stochastic parameter perturbation of this study, MTSP, with a Markov process time correlation and also a spatial correlation of nodes with common neighbors, shows a strong impact. The difference of the ensemble mean compared to the reference run, as well as the generated ensemble spread, is of considerable magnitude. Both increase during long-lasting periods of convergent driftand constant ice coverage. This reflects the cumulative nature of the deformation processes in the sea ice cover. On the other hand, the influences of melting and freezing as well as divergent driftpatterns lead to a decrease in the ensemble spread and also in the mean influence of the P* perturbations. Because of the lack of appropriate data, especially in regions with multiyear ice where the influence of the perturbation schemes is strongest, we cannot derive an exact calibration of spread compared to the ensemble mean error, either for sea ice thickness or for concentration. Compared to the model's interannual spread in sea ice thickness, though, the ensemble spread is of comparable magnitude. Using this finding as a first indication, the ensemble spread seems to be reasonable. The spread of sea ice thickness within the ensemble produced by the various P* perturbations applied may indicate a range of possible P* values for future data assimilation studies. It shows the general influence local changes of P* can have on the sea ice thickness distribution patterns within the model. Within this study an appropriate range of values seems to be P* 2 (5000Nm22, 35 000Nm22). To summarize the results of this paper, first of all it would be good to have a more extensive dataset for sea ice thickness, to be able to evaluate (not only) stochastic sea ice parameterization schemes more thoroughly. Long-time measurements at fixed locations such as ULS data are of great value for this purpose but are hard to obtain in regions with a thick perennial ice cover. Nevertheless, even without such observations our results prove that it is of importance to look at the inherent uncertainties of parameterizations. One approach to do so is to include stochasticity in the formulation of parameterizations, as Palmer (2012) suggests in view of a probabilistic Earth-system simulator. The effect of parameterization uncertainties is going to be even stronger in fully coupled atmosphere-ocean-sea ice models and will influence the potential of the model to predict future states of the climate system. We are currently working on the inclusion of stochastic sea ice parameterizations in such a fully coupled model. Even though we have as yet only presented a first approach to stochastic perturbations in a sea ice model, we think that such probabilistic parameterizations should be implemented in future climate ensemble prediction systems to include the model uncertainty at a level where it first appears in the equations. This is especially of importance when it comes to seasonal and interannual predictions. In the process of developing parameterizations and assessing their inherent uncertainties it is essential to think about model resolution, as increased model resolution has an impact on the statistics the parameterizations are based upon. These changing statistics might be more readily applied to a stochastic parameterization than a purely deterministic one. Perturbing P* in a sea ice model is a first step toward addressing these topics and serves as a tool to test and evaluate the sensitivities of the model and the mesh it uses. Even though a constant P* value might be sufficient for many problems, it should be kept in mind what an influence slightly perturbed parameters may have, even if the perturbations are symmetric. Further work will now be directed toward other perturbation schemes, especially concerning spatial correlation, and toward other parameters of the sea ice model, especially toward thermodynamic parameters such as the sea ice albedo. In a first step, we compared our approach of a stochastically perturbed P* with a realization of the stochastically perturbed parameterization tendencies (SPPT) scheme by Buizza et al. (1999). Here, the parameterization of the internal forces Fint in Eq. (1) has been perturbed by multiplicative noise. In our experiment with the SPPT scheme there was no spatial correlation between nodes and the reassignment time steps used were similar to those of our reassignment time step perturbation approach, RP. The results of this first application of the SPPT scheme to sea ice dynamics showed changes of similar magnitude as theRP approach (not shown). Therefore implementation of temporal correlation in the Markov process time correlated approach, MTP, and implementation of temporal and spatial correlation in the MTSP approach proved to be much more influential concerning sea ice dynamics than this specific realization of the SPPT scheme. We have also started to apply the general idea of perturbing traditionally fixed parameters in some preliminary experiments to the parameter C in Eq. (7) of the sea ice rheology and also to the lead closing parameter h0 used in the parameterizations of the continuity equation for the sea ice concentration. The latter parameter influences the relation between increase in sea ice thickness and increase in sea ice concentration (lead closing) during freezing conditions. First results have shown that parameter perturbations concerning C are of less importance compared to P*, whereas the perturbations of h0 show quite a strong response with respect to the sea ice concentration in the Southern Hemisphere (not shown). For Southern Ocean sea ice, the influence of P* perturbations is small, which is partly due to the small amount of modeled multiyear ice in this FESOM version, but probably also reflects the natural differences between Arctic and Antarctic. Further studies connected to these first results will be directed toward the different influences parameter perturbations have on changes in ice thickness and concentration distribution in the Northern and Southern Hemisphere, which have shown quite diverging results so far. Another issue that will be investigated in the future is the influence of methods that can be used to include spatial correlation of random values. It is of great interest to see how this correlation, and therefore the entire stochastic parameterization, is influenced by the varying mesh resolution in a finite element model. For investigations of the resolution-dependent behavior of stochastic and also deterministic parameterizations FESOM is an excellent tool. Acknowledgments. Wethank SergeyDanilov, Thomas Jung, and Martin Losch for their help and advice in preparing this publication. Our thanks go to Hannelore Witte for providing the Greenland Sea and Fram Strait ULS data and to Stefan Hendricks, who helped us with his knowledge of Arctic sea ice thickness distributions. Finally, we thank the two anonymous reviewers for their helpful comments and advice. The Beaufort Sea ULS data and the sea ice concentration data were obtained from the National Snow and Ice Data Center (NSIDC), Boulder, Colorado (http://www.nsidc.org/). The sea ice thickness ICESat-data were obtained fromtheNASAJet Propulsion Laboratory (http://rkwok.jpl.nasa.gov/). REFERENCES Bardsley, W. E., 2006: Note on y-truncation: A simple approach to generating bounded distributions for environmental applications. Adv. Water Resour., 30, 113-117. Berner, J., G. J. Shutts, M. Leutbecher, and T. N. Palmer, 2009: A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF Ensemble Prediction System. J. Atmos. Sci., 66, 603-626. Bright, D. R., and S. L. Mullen, 2002: Short-range ensemble forecasts of precipitation during the Southwest monsoon. Wea. Forecasting, 17, 1080-1100. Buizza, R., M. Miller, and T. N. Palmer, 1999: Stochastic representation of model uncertainties in the ECMWF Ensemble Prediction System. Quart. J. Roy. Meteor. Soc., 125, 2887-2908. Cavalieri, D., C. Parkinson, P. Gloersen, and H. J. Zwally, 1996: Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I passive microwave data (updated yearly).National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice. gd.html.] Flato,G.M., 1998: The thickness variable in sea-icemodels. Atmos.- Ocean, 36, 29-36. Harder, M., 1996: Dynamics, roughness, and age of Arctic sea ice- Numerical investigations with a large-scale model. Alfred Wegener Institute for Polar andMarine Research Rep. 203, 127 pp. _____, and H. Fischer, 1999: Sea ice dynamics in the Weddell Sea simulated with an optimized model. J. Geophys. Res., 104 (C5), 11 151-11 162. Hibler, W. D., 1979: A dynamic thermodynamic sea ice model. J. Phys. Oceanogr., 9, 815-846. Hunke, E. C., and J. K. Dukowicz, 1997: An elastic-viscous-plastic model for sea ice dynamics. J. Phys. Oceanogr., 27, 1849-1867. _____, and W. H. Lipscomb, 2010: CICE: The Los Alamos Sea Ice Model: Documentation and software user's manual, version 4.1. Tech. Rep. LA-CC-06-012, T-3 Fluid Dynamics Group, Los Alamos National Laboratory, 76 pp. Kloeden, P. E., and E. Platen, 1992: Numerical Solution of Stochastic Differential Equations. Springer-Verlag, 636 pp. Kotecha, J. H., and P. M. Djuric, 1999: Gibbs sampling approach for generation of truncated multivariate Gaussian random variables. Proc. ICASSP 1999, Phoenix, AZ, IEEE, Vol. 3, 1757-1760. Kwok, R., G. F. Cunningham, H. J. Zwally, and D. Yi, 2007: Ice, Cloud, and Land Elevation Satellite (ICESat) over Arctic sea ice: Retrieval of freeboard. J. Geophys. Res., 112, C12013, doi:10.1029/2006JC003978. _____, _____, M. Wensnahan, I. Rigor, H. J. Zwally, and D. Yi, 2009: Thinning and volume loss of the Arctic Ocean sea ice cover: 2003-2008. J. Geophys. Res., 114, C07005, doi:10.1029/ 2009JC005312. Large, W. G., and S. G. Yeager, 2009: The global climatology of an interannually varying air-sea flux data set. Climate Dyn., 33, 341-364. Lemke, P., E. W. Trinkl, and K. Hasselmann, 1980: Stochastic dynamic analysis of polar sea ice variability. J. Phys. Oceanogr., 10, 2100-2120. Li, X., M. Charron, L. Spacek, and G. Candille, 2008: A regional ensemble prediction system based on moist targeted singular vectors and stochastic parameter perturbations. Mon. Wea. Rev., 136, 443-462. Lin, J. W.-B., and J. D. Neelin, 2000: Influence of a stochastic moist convective parameterization on tropical climate variability. Geophys. Res. Lett., 27, 3691-3694. Melling, H., and D. A. Riedel, 2008: Ice draftand ice velocity data in the Beaufort Sea, 1990-2003. National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http:// nsidc.org/data/g02177.html.] Owens, W. B., and P. Lemke, 1990: Sensitivity studies with a sea ice-mixed layer-pycnocline model in the Weddell Sea. J. Geophys. Res., 95 (C6), 9527-9538. Palmer, T. N., 2012: Towards the probabilistic Earth-system simulator: A vision for the future of climate and weather prediction. Quart. J. Roy. Meteor. Soc., 138, 841-861. Sidorenko, D., Q. Wang, S. Danilov, and J. Schröter, 2011: FESOM under coordinated ocean-ice reference experiment forcing. Ocean Dyn., 61, 881-890. Stössel, A., P. Lemke, and W. B. Owens, 1990: Coupled sea icemixed layer simulations for the Southern Ocean. J. Geophys. Res., 95 (C6), 9539-9555. Timmermann, R., S. Danilov, J. Schröter, C. Böning, D. Sidorenko, and K. Rollenhagen, 2009: Ocean circulation and sea ice distribution in a finite element global sea ice-ocean model. Ocean Modell., 27, 114-129. Tremblay, L.-B., and M. Hakakian, 2006: Estimating the sea ice compressive strength from satellite-derived sea ice driftand NCEP reanalysis data. J. Phys. Oceanogr., 36, 2165-2172. Wang, Q., S. Danilov, and J. Schröter, 2008: Finite element ocean circulation model based on triangular prismatic elements, with application in studying the effect of topography representation. J. Geophys. Res., 113, C05015, doi:10.1029/2007JC004482. Witte, H., and E. Fahrbach, 2005: AWI moored ULS data, Greenland Sea and Fram Strait, 1991-2002. National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http://nsidc.org/data/g02139.html.] Yu, Y., G. A. Maykut, and D. A. Rothrock, 2004: Changes in the thickness distribution of Arctic sea ice between 1958-1970 and 1993-1997. J. Geophys. Res., 109, C08004, doi:10.1029/ 2003JC001982. STEPHAN JURICKE, PETER LEMKE, RALPH TIMMERMANN, AND THOMAS RACKOW Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany (Manuscript received 28 June 2012, in final form 19 November 2012) Corresponding author address: Stephan Juricke, Alfred Wegener Institute, Bussestraße 24, F-406, 27570 Bremerhaven, Germany. E-mail: [email protected] APPENDIX A Bounded Distributions For this study bounded distribution functions for the random numbers of the perturbations are needed. In Fig. A1 the corresponding distribution functions of the so-called x- and y-truncated Gaussian distributions (Bardsley 2006) can be seen, with a fixed variance s2 and zero mean of the initial Gaussian distribution. An additional bounded, Gaussian-like random number can be generated by the transformation ... (A1) where y(i, j) is a random number from a Gaussian distribution, again with variance s2 and zero mean, and a is the bounding parameter. In this transformation b is some fixed value that meets the inequality ... with b 5 1.4/slim used here. The value b is necessary to keep the maximum of the transformed distribution function at 0. The quantity slim is a value greater or equal to the standard deviation s of the Gaussian distribution that is going to be transformed. The distribution function after the transformation of a Gaussian distributed random number y(i, j) is also shown in Fig. A1. APPENDIX B Derivation of the Markov Process Time Correlated Perturbation Equation (9) in section 3b can be derived by discretizing a stochastic differential equation of the form ... (B1) using the Euler-Maruyama scheme (Kloeden and Platen 1992), where (Wt)t$0 is a Wiener process and s2 W 5s2/(Dt)pert. Under these conditions the discretization of Eq. (B1) yields ... (B2) with ... where t is the relaxation time, as already mentioned in section 3b. The mean of a random variable y at node i described by Eq. (9) is zero for all time steps j, while the variance is given by ... (B3) and depends on the amount of time steps the node has been covered by ice since it was last ice-free, ( j2 j0), with j being the current time step and j0 the time step when the node was last ice-free. The resulting random variable is Gaussian distributed because all the summands are Gaussian distributed. For ... the variance is limited by ... (B4) For the transformation of the random number y(i, j) into the bounded x(i, j) used in Eq. (9) the simple Eq. (A1) can be used, or corresponding transformations for the cases of x- and y-truncated Gaussian distributions. Kotecha and Djuric (1999) gave such a transformation for the x-truncated Gaussian distribution. In case of the more complicated y-truncated Gaussian distribution the corresponding transformation needs to be implemented via, for instance, Newton's method. For both transformations the variance at the current time step is needed. For the simple transformation via Eq. (A1) the value slim is set to the limiting variance in Eq. (B4) for all time steps. (c) 2013 American Meteorological Society [ Back To TMCnet.com's Homepage ]