In this paper a method to directly assimilate brightness temperatures from the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E) to produce ice concentration analyses within a three-dimensional variational data assimilation system is investigated. To assimilate the brightness temperatures a simple radiative transfer model is used as the forward model that maps the state vector to the observation space. This allows brightness temperatures to be modeled for all channels as a function of the total ice concentration, surface wind speed, sea surface temperature, ice temperature, vertically integrated water vapor, and vertically integrated cloud liquid water. The brightness temperatures estimated by the radiative transfer model are sensitive to the specified values for the sea ice emissivity. In this paper, two methods of specifying the sea ice emissivity are compared. The first uses a constant value for each polarization and frequency, while the second uses a simple emissivity parameterization. The emissivity parameterization is found to significantly improve the fit to the observations, reducing both the bias and the standard deviation. Results from the assimilation of brightness temperatures are compared with those from assimilating a retrieved ice concentration in the context of initializing a coupled ice–ocean model for an area along the east coast of Canada. It is found that with the emissivity parameterization the assimilation of brightness temperatures produces ice concentration analyses that are in slightly better agreement with operational ice charts than when assimilating an ice concentration retrieval, with the most significant improvements during the melt season.
An accurate estimate of the sea ice state is critical for providing information to ensure safe ship navigation in ice-infested waters, for improved numerical weather prediction (NWP) near ice-covered regions, and for climate studies. If the Arctic continues to warm as projected, there will be an increased need for accurate sea ice information because of an increase in ship traffic for transport and natural resource extraction in ice-covered regions. The recent results of the Concordiasi project demonstrate the benefit of assimilating satellite data over sea ice (Bouchard et al. 2010) to NWP, while the positive impact an improved sea ice state can have on weather forecasts is evident in the coupled Gulf of St. Lawrence prediction system, now running operationally at the Canadian Meteorological Centre (Pellerin et al. 2004).
Passive microwave sensors are useful for monitoring the state of sea ice because of their ability to operate in cloudy conditions and without sunlight. This has led to the development of numerous methods to estimate ice concentration from passive microwave brightness temperatures. A comprehensive review of these methods can be found in Andersen et al. (2006). Ice concentration derived from passive microwave data is a frequent choice of observation type used in sea ice data assimilation (Stark et al. 2008; Lindsay and Zhang 2006; Lisaeter et al. 2003). Lisaeter et al. (2003) assimilate ice concentration retrieved from the Special Sensor Microwave Imager (SSM/I) using an ensemble Kalman filter (EnKF). They found that their method improved estimates of ice concentration relative to a free-run coupled ice–ocean model experiment, but that the background error variance was underpredicted, which they attributed to bias in the model that was not accounted for. They demonstrated the complexity of the spatial and multivariate covariances of the background error. Stark et al. (2008) assimilated ice concentration and ice motion derived from passive microwave sensors. They used an optimal interpolation scheme and a novel method to incorporate the ice velocity increment through adding a stress term to the ice momentum equation. This was necessary to keep the model from rejecting the velocity increments due to the short time scales associated with the elastic–viscous–plastic ice rheology. They carried out their assimilation over an entire ice season for both the Arctic and Antarctica and found that assimilating data led to significant improvements, especially during the summer. They had some difficulty in getting the model to retain thin ice increments, possibly because there were no updates to the ocean temperature. In contrast, the study by Caya et al. (2010) used a three-dimensional variational data assimilation system (3D-Var) with multivariate background error covariances that allowed the ocean temperature and salinity to be updated when the ice concentration is changed by the assimilation. However, they did not find that including the update to the oceanic variables had a significant impact on their results, possibly because the background error was modeled as spatially and temporally homogeneous. With such an approach, it is not possible to account for the complexity of the background error covariances, such as that found in the study by Lisaeter et al. (2003).
The approach used here differs from others presented previously in that instead of assimilating the ice concentration retrieved from passive microwave brightness temperatures, the observed brightness temperatures are assimilated directly. In the NWP community it has been found that directly assimilating satellite brightness temperatures, rather than retrieved temperature and moisture profiles, has a significant positive impact on weather forecasts (Saunders et al. 1999; Derber and Wu 1998). However, the work presented here differs from past work in the NWP community in that the low-frequency channels, which are sensitive to the surface conditions over the ocean where sea ice may be present, are being assimilated. The use of surface sensitive observations other than those over the ice-free ocean has only recently been a topic of research in NWP, with work being led by Karbou et al. (2005, 2006) for the assimilation of surface-sensitive channels from the Advanced Microwave Sounding Unit (AMSU) over land. They compared using the operational system, which assigns an emissivity value based on surface type, to using a retrieved emissivity (i.e., one calculated with a radiative transfer equation using passive microwave observations). They calculated a retrieved emissivity at an atmospheric window channel, and then assigned this emissivity to the sounding channels. They noted a reduction in the predicted minus observed (P − O) brightness temperature bias and standard deviation when using the retrieved emissivity, with a more significant improvement when they used emissivity calculated each day as compared to using an average over the previous month. A very similar method was also used by Bouchard et al. (2010) to assimilate AMSU observations over Antarctica. The method used by Bouchard et al. (2010) calculates the emissivity of the surface without regard as to whether the surface is ice or open water. This is because their primary interest is in assimilating surface-sensitive channels to assist in estimating the atmospheric state.
The goal of the present study is to estimate surface conditions, in particular the sea ice concentration, by directly assimilating passive microwave brightness temperatures. The emissivity of sea ice is calculated using a radiative transfer model and a background estimate of the ice concentration, ice temperature, and other variables that affect the brightness temperatures. The emissivity is calculated using only a subset of the data where the ice concentration is high and atmospheric effects are small. The calculated emissivity is correlated with the ice thickness for thin ice, and an average emissivity is calculated for thicker ice. The emissivities are calculated for each frequency, polarization, and for each month. Results from experiments that directly assimilate brightness temperatures using either the calculated emissivity or fixed values of emissivity for first-year ice (Eppler et al. 1992) interpolated to the frequencies of the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E), will be compared together with the more conventional method of assimilating a retrieved ice concentration. The assimilation system and method used to calculate the emissivity are described in section 2. This is followed by a demonstration of the impact of the emissivity parameterization in section 3, and then a comparison between the direct assimilation of brightness temperatures and assimilation of a retrieved ice concentration in section 4. A verification of the sea ice concentration and thickness is given in section 5, followed by the conclusions in section 6.
a. Ice–ocean model
The ice–ocean model used in this study is the Community Ice–Ocean Model (CIOM) that covers the region off the east coast of Canada, described by Yao et al. (2000). It consists of a multicategory ice model with the vertical sigma coordinate system Princeton Ocean Model (Blumberg and Mellor 1987). The model resolution is in latitude and in longitude. Prognostic equations are solved for the three-dimensional ocean velocity, potential temperature, salinity, and turbulence quantities on 16 sigma levels. The ice model describes both ice dynamics and thermodynamics, and consists of a diagnostic equation for the ice velocity, an equation governing the heat transfer through the ice, prognostic equations for 24 ice-thickness categories, and uses an elastic–plastic rheology. The current version of the model differs from that described in Yao et al. (2000) in the following ways: the shortwave radiation is taken from Shine et al. (1984) with the cloud factor from Laevastu (1960) following the study by Dunlap et al. (2007), the ocean boundary conditions evolve continuously from monthly means, the ice thickness at the northern boundary evolves based on local thermodynamics, and the model can be restarted precisely from a previous run. Atmospheric forcing to CIOM is provided by the Canadian Global Environmental Multiscale (GEM) model (Côté et al. 1998) short-term forecasts, and are refreshed every 3 h. CIOM has been found to produce ocean currents and ice thickness in general agreement with known values.
b. Assimilation system
The method used to assimilate the observations is a 3D-Var method that minimizes the difference between the analysis and the background, and the analysis and the observations. The standard 3D-Var cost function has the following form:
where xb is the background state (a 12-h model forecast for the ice and ocean analysis variables), y is the vector of observations, H is the (possibly nonlinear) observation operator that maps the state vector to the observation space, is the background error covariance matrix, and is the observation error covariance matrix. The state x that minimizes the cost function is referred to as the analysis. For practical reasons, the cost function is transformed into its incremental form (Courtier et al. 1994) and preconditioned with respect to the background error term. This results in the cost function:
where ξ is the control vector and is the observation operator linearized about xb. In deriving (2) the approximation has been made. This approximation results in J being a quadratic function with respect to ξ. The analysis is related to the control vector by
c. Assimilation of retrieved ice concentration
The first type of observation assimilated is the ice concentration calculated from brightness temperatures measured by the AMSR-E sensor using the National Aeronautics and Space Administration (NASA) Team 2 (NT2) algorithm (Markus and Cavalieri 2000). This algorithm uses a set of precomputed tie points for first-year ice, multiyear ice, thin ice, and layered ice, and the polarization and gradient ratios between the brightness temperatures from the horizontally and vertically polarized 18.7-, 36.5-, and 89.0-GHz channels. Because the algorithm uses ratios of brightness temperatures, it is less sensitive to ice temperature than algorithms that use the brightness temperatures directly. Known problems with NT2 include: weather contamination, for example mistaking the wind-roughened sea surface for ice; biases caused by the presence of meltponds, where meltwater cannot be distinguished from open water; and, mistaking thin ice for open water, due to the fact that the emissivity of thin ice is close to that of open water. NT2 ice concentration retrievals using AMSR-E data have been validated against those from Landsat imagery (Cavalieri et al. 2006). Reported differences ranged between 1% and 16% over consolidated first-year ice and over a diffuse marginal ice zone, respectively. Previous studies of the NT2 algorithm using SSM/I data estimated the error to be between 5% and 17% (Markus and Dokken 2002; Shokr and Markus 2006). In the current study, the ice edge is very diffuse in the southern portion of the domain and an observation error standard deviation of 20% to specify for the retrieved NT2 ice concentration is therefore used. When the retrieved ice concentration is assimilated, the observation operator H is a linear interpolation operator that interpolates the background and analysis total ice concentration (computed from the partial concentration for all thickness categories) to the observation times and locations.
The background error covariances for the partial ice concentrations for each of the models ice thickness categories, ocean temperature, and ocean salinity were determined using an ensemble Kalman filter, as described in Caya et al. (2010). To avoid the creation of spurious ice seen in preliminary experiments, the ocean variables in the model are not updated by the 3D-Var analysis in the current study. Consequently, only the ice variables in the model are modified when initializing the 12-h forecast used to produce the background state for the next analysis.
d. Assimilation of brightness temperatures
The second type of observations assimilated are the AMSR-E brightness temperatures, by using a radiative transfer model as part of the observation operator in (2). The brightness temperatures at five frequencies, 6.9, 10.7, 18.7, 23.8, and 36.5 GHz, and both horizontal and vertical polarizations are assimilated. The 89-GHz channels are not being assimilated because they are more sensitive to the atmospheric conditions than the lower-frequency channels. In the study by Andersen et al. (2006), it was found that this sensitivity reduced the accuracy of retrieval algorithms when the 89-GHz channel was used in regions close to the ice edge. Since the eastern ice edge is a particularly dominant feature in our domain, it was decided to exclude the 89-GHz channel for this study. It is expected that the 89-GHz channels would provide useful information at colder temperatures and away from the ice edge where there are less atmospheric effects, for example, during the winter in the Arctic. This will be addressed in a future study.
The radiative transfer model (RTM) developed by Wentz and Meissner (2000) is used as the starting point for the calculation of brightness temperatures. This model calculates the upwelling brightness temperature over the open ocean at polarization p and frequency ν, seen by the satellite according to
where Tu is the upwelling atmospheric emission, TD is the downwelling atmospheric emission, Tow is the sea surface temperature, εow is the sea surface emissivity, τ is the atmospheric transmissivity, Ωow is the open-water surface scattering term, and TC is the cold space temperature. The upwelling and downwelling atmospheric emission, sea surface emissivity, scattering term, and atmospheric transmissivity are calculated within the RTM as a function of the 10-m wind speed, the integrated water vapor, the integrated cloud liquid water, and the surface temperature. The background values of the atmospheric variables are obtained from short-term GEM forecasts. The background value of the sea surface temperature is taken from the ice–ocean model.
To accommodate sea ice, the calculation of the upwelling brightness temperature is modified to be
where Cice, Tice, and εice are the ice concentration, ice temperature, and ice emissivity, respectively. The background values of ice concentration and ice temperature are obtained from the ice–ocean model; the sea ice emissivity must be specified.
Unlike open water, where emissivity models exist and provide reliable estimates of the surface emissivity, similar models for sea ice are still under development (Tonboe 2010; Guedj et al. 2010). Ice emissivity is a complicated function of the surface conditions; snow cover and layering, surface moisture, ice brine volume, and surface roughness, among others (Stroeve et al. 2006; Tonboe 2010; Guedj et al. 2010; Haggerty and Curry 2001; Comiso et al. 1989; Shokr et al. 2009). Without detailed knowledge of the surface conditions it is difficult to accurately model the emissivity. Here a detailed model of the sea ice emissivity is not attempted, but instead the radiative transfer model is used to calculate the emissivity using a subset of the observations. A simple parameterization of the ice emissivity in terms of the ice thickness is developed and used for one of the experiments. This approach was motivated by the study of Naoki et al. (2008) where it was shown that an approximately linear relationship exists between ice thickness and emissivity for thin ice (ice of thickness less than ≈0.3 m). It was preferable to start with a simple linear fit and determine if the emissivity parameterization had an impact on the assimilation, rather than spend more time refining the emissivity parameterization. This could be addressed in a future study. In Naoki et al. (2008) it was found that the sensitivity of the emissivity to sea ice thickness is greater for the low frequencies than for the high frequencies, and that the variation in emissivity due to polarization is greater than that due to the frequency. This has also been found in the present study. The correlation between ice thickness and emissivity for the different frequencies and polarizations are given in Table 1. It should be noted that since multiyear ice is not prevalent within the domain of interest, it is only necessary to consider the emissivity value for first-year ice.
To calculate the ice emissivity, the same equation that is used to calculate the brightness temperature [(5)] is used, but instead of solving for brightness temperature the observed brightness temperature is used and the equation is solved for the sea ice emissivity. This calculation is performed separately for each month using data from every fourth day. Only points where the atmospheric transmissivity is greater than 0.8 for the frequency for which the calculation is being carried out are included to minimize errors. In addition, during winter months only points where the ice concentration is greater than 0.9 are used, while during the spring it was necessary to decrease this threshold to 0.8 for April and 0.7 for May and June, because of the reduced number of points in the domain with high ice concentration. The data is split into two categories, depending on whether the ice thickness is greater or less than 0.3 m. For the data where the thickness is less than 0.3 m, a least squares fit of the emissivity using a linear function of ice thickness is carried out. For the data where the thickness is greater than 0.3 m, the average emissivity is calculated. For these calculations, the ice thickness is taken from the ice model in the experiment that assimilates the NT2 ice concentration retrieval. As an example, a scatterplot of the emissivity versus ice thickness calculated for the horizontally polarized 6.9-GHz channel during the month of December is shown in Fig. 1. At this time almost the entire domain is composed of thin ice, and the correlation between the ice thickness and emissivity is strong. Correlation coefficients between ice thickness and emissivity (for thickness less than 0.3 m) are given in Table 1 for all channels. After March, there is very little thin ice in the domain, and therefore only the average ice emissivities are calculated for these months. It can be seen that the correlation is higher for the low-frequency channels, and also for the horizontally polarized channels. The values of emissivity calculated as a function of ice thickness were found to be in good agreement with those given in Naoki et al. (2008), which were determined during a field campaign where ice thickness, brightness temperatures, and ice temperature were measured simultaneously.
The average emissivities for the thick ice (December–March) and for all ice (April–June) for each month are shown in Fig. 2. The computed values are compared against the values from Eppler et al. (1992) for first-year ice because these are the values that are used in the development of the NT2 algorithm. It can be seen that the emissivities calculated using the RTM and those from Eppler et al. (1992) are in reasonably good agreement, considering that the radiative transfer model does not take into account many of the factors affecting emissivity at these frequencies, such as surface roughness. Surface roughness tends to lower the emissivity in the horizontally polarized channels, and increase it very slightly in the vertically polarized channels Stroeve et al. (2006). It should also be noted that the emissivity values in Eppler et al. (1992) are computed using data from surface and aircraft-based sensors, while our values are calculated using satellite data that have a much larger footprint. For the satellite data, small openings in the ice containing open water would not be resolved, and would act to lower emissivity of the sea ice because the emissivity of open water is much lower than that of ice. This would have a larger impact for the horizontally polarized channels than the vertically polarized channel because the open-water emissivity is lower for the horizontally polarized channels.
In this study the two lowest-frequency channels, 6.9 and 10.7 GHz, being assimilated are normally not used in ice concentration retrieval algorithms because of their coarse resolution. However, these channels are beneficial to the analysis because of their high sensitivity to sea ice, and low sensitivity to atmospheric effects. The sensitivity of the assimilated AMSR-E channels to the atmosphere and to the ice concentration is shown in Fig. 3, which shows the component of the normalized Jacobian operator, , which corresponds with the sensitivity to ice concentration. This is shown both for a case with very little atmospheric effects (water vapor of 3 mm and cloud liquid water of 0.002 mm), and for a case with significant atmospheric effects (water vapor of 20 mm and cloud liquid water of 0.8 mm). In both cases the low-frequency channels are more sensitive to the ice concentration than the high-frequency channels, and the sensitivity of the higher-frequency channels to ice concentration is reduced in the presence of significant atmospheric effects.
The footprint size of the different channels used ranges from 43 km × 75 km for the 6.9-GHz channels, to 14 km × 8 km for the 36.5-GHz channels. To reduce land contamination, any observation within a distance of half of the maximum dimension of the elliptical footprint from the land boundary was not assimilated. Because the footprint sizes are different for each frequency, this criterion was applied to each frequency separately. Consequently, only the higher-frequency channels are assimilated close to land. For the two lowest-frequencies assimilated (6.9 and 10.7 GHz) the average radius of the sensor footprint is larger than the grid spacing of the 3D-Var analysis and background states. For these frequencies the background and analysis points within the sensor radius are averaged as part of the observation operator to give the value at the observation location. A top-hat filter is used with radius equal to the radius of a circle with the same area as the area of the elliptical footprint. A different radius is used for each frequency. With this methodology the low-frequency channels will only influence the average analysis increment over their respective footprints, while the higher-frequency channels will be free to influence the smaller-scale component of the analysis increment. This is in contrast to the NT2 algorithm that combines frequencies having different sizes of footprints to calculate a single estimate of the ice concentration, which can lead to smearing of the ice edge (Hwang and Barber 2006).
To be able to use the radiative transfer model the analysis vector includes the ice concentration in each ice thickness category (from which the total ice concentration is computed), ocean surface temperature, surface wind speed, vertically integrated water vapor, vertically integrated cloud liquid water, and ice temperature taken at a depth of 1.8 cm. The depth of 1.8 cm was chosen because it represents a nominal penetration depth at 17 GHz (Mathew 2007). The background error covariances for the partial ice concentrations for each of the models ice thickness categories, and ocean temperature are the same as those used in the assimilation of the retrieved ice concentration. To estimate the background error standard deviations for the three atmospheric variables, the standard deviation of the difference between the background state from the GEM forecast fields, and the same variables estimated by AMSR-E retrievals (Wentz and Meissner 2000), were calculated over the ice season. These values, which most likely overestimate the real uncertainty in the GEM fields, are given in Table 2. For ice temperature, the background error standard deviation was arbitrarily chosen to be 5 K for all experiments reported here. The sensitivity of the results to the ice temperature variance was checked by repeating the assimilation of brightness temperatures with the ice temperature background error variance doubled and halved. The resulting ice concentration analyses and forecasts were found to be insensitive to these changes to the ice temperature background error variance. For the variables given in Table 2 the background errors were assumed to be uncorrelated with all other analysis variables. Though this assumption is not likely accurate, it is a reasonable starting point. Note that although the analysis increments for the ocean temperature, ice temperature, and atmospheric variables are not applied to correct the background state when initializing the subsequent short-term model forecast, it is important that they are included as analysis variables. For example, by including the atmospheric variables they can, when appropriate, explain the part of the innovation, which is the difference between the observed brightness temperatures y, and the brightness temperatures calculated using the background state H(xb), which is due to atmospheric effects. This reduces the appearance of spurious ice over the open ocean that arises because of the difference between the observed brightness temperatures y, and those calculated using the background state H(xb).
The observation error covariance matrix must also be specified for the assimilated brightness temperature channels. To estimate , the variance of the innovation is used. This method of calculating the observation error represents an approximate upper bound on the observation error, assuming that the observation and background errors are uncorrelated (Desroziers et al. 2005). The upper bound was used partly because the matrix was assumed diagonal, although there is in reality some interchannel and spatial correlation due to both instrument errors and especially forward model errors, including errors in the emissivity. To calculate the innovation, the background values of ice concentration and ice temperature were taken from the experiment in which ice concentration retrievals were assimilated. The variance for each channel was calculated separately for each month.
e. Experimental setup
Results from three experiments are presented in this study. In the first experiment the ice concentration retrieval calculated using AMSR-E brightness temperatures and the NT2 algorithm are assimilated as described in section 2c. The second and third experiments assimilate brightness temperatures with the constant sea ice emissivity from Eppler et al. (1992), and with our monthly-varying sea ice emissivity parameterization, respectively.
The ice concentration, sea surface temperature, and ice temperature, produced by the first experiment that assimilates the NT2 ice concentration, along with the integrated water vapor, integrated cloud liquid water, and wind speed from the GEM forecasts, are used to calculate predicted brightness temperatures over the entire ice season using two separate versions of the radiative transfer model. The first version uses sea ice emissivity values from Eppler et al. (1992) in the radiative transfer model, while the second uses our emissivity parameterization. There are two purposes for calculating these two sets of brightness temperatures. The first is to see if the emissivity parameterization improves the fit between the predicted and observed brightness temperatures, and the second is to compute the innovation statistics. The innovation statistics are then used to calculate the observation error covariance matrix used to assimilate brightness temperatures in experiments 2 and 3.
The brightness temperatures and retrieved ice concentration are assimilated every 12 h, with the analysis time at the middle of the 12-h window, at 0600 and 1800 UTC. The GEM output, which is available at 3-h intervals, is linearly interpolated in time and space to the AMSR-E observation time and location. Because there is a spinup time needed to develop realistic clouds in the model, only forecasts with at least a 6-h lead time were used. The assimilation is carried out over the same time period used by Caya et al. (2010): 5 December 2006–30 June 2007. For all experiments, the quality control used here rejects observations for which no ice occurs in the background state within 100 km. The domain chosen is the east coast of Canada (Fig. 4). The verification data chosen are daily ice charts for the Newfoundland region produced by the Canadian Ice Service (Carrieres et al. 1996). These charts are produced manually by analysts who subjectively combine recent history, climatology, and model forecast with observations from various sources (e.g., synthetic aperture radar, visual and infrared sensors, passive microwave imagery, as well as ship and aircraft reports). The relative weight given to passive microwave observation is small, and hence these analyses can be considered to be independent of the data assimilated here. The dark shaded region in Fig. 4 is the verification domain, chosen to be the weekly maximum ice extent observed over the past 30 years, south of 57°N.
3. Impact of the emissivity parameterization on the ability to fit observations
In this section the statistics of the P − O brightness temperatures calculated using the radiative transfer model with the background state from the experiment that assimilated the retrieved ice concentration (experiment 1) are compared. Two sets of results are shown, one with the constant emissivity, and one with the emissivity parameterization, for which the parameters are computed separately for each month. This will tell us to what extent the emissivity parameterization improves the fit to the observations.
In the development of the 3D-Var system it is assumed that the background and observation errors are unbiased and follow a normal distribution. However, even if this is the case, nonlinearity of the observation operator (i.e., the radiative transfer model) will cause the innovations to be nonnormal and possibly biased. In addition, because of possible biases in the background state of the surface temperature and emissivity, the biases in the innovation for microwave channels that are sensitive to the surface can be significant (Karbou et al. 2006; Bouchard et al. 2010), which can lead to biased analysis increments. It is expected that the emissivity parameterization scheme should reduce the bias because it is calculated using a subset of the same background and observation values used to compute the innovation statistics. For example, to some extent biases in the surface temperature will be compensated for in the calculated emissivity. As an example, a histogram of the innovations over the month of December for the 6.9-GHz horizontally polarized channel is shown in Fig. 5, where it can be seen that using the emissivity parameterization reduces the bias and the large positive tail, and brings the innovations closer to a normal distribution. The bias, standard deviation, and skew for this example are 14.9, 24.3, and 2.28 K, respectively, with the constant emissivity, and 1.8, 12.4, and 1.84 K with the emissivity parameterization. The excess kurtosis, which is a measure of how large the tails of the distribution are relative to a normal distribution, is 7.07 with the constant emissivity and 4.57 for the emissivity parameterization, with the lower value of excess kurtosis indicating a distribution that is closer to normal.
In Fig. 6 the monthly averages of the bias and standard deviations for all channels are plotted. It can be seen that when using the constant emissivity, the bias, indicated by the thick lines, for the horizontally polarized channels is very high in December (11–17 K) and again during May and June (4–9 K). This is probably due to the presence of thin ice in December, and meltwater on the surface during May and June, both of which are known to reduce the surface emissivity. With the emissivity parameterization the bias is significantly reduced. During the remaining months (January, February, March, and April) the impact of the emissivity parameterization is neutral, with the negative bias being −2 to −8 K whereas there is a positive bias of the same magnitude when a constant emissivity is used. For the vertically polarized channels the impact of the emissivity parameterization on the bias is less pronounced. The positive bias in December, May, and June is reduced, as is the negative bias during the months of January, February, March, and April.
The monthly averages of the standard deviation are shown as the dotted lines in Fig. 6. For the vertical channels it can be seen that with the emissivity parameterization the standard deviation is fairly constant across the season, and is around 10 K. With the constant emissivity the standard deviation is higher during December and January, being between 10 and 15 K. For the horizontal channels the standard deviation during December and January is significantly reduced (from 20–30 to 10–15 K) by using the emissivity parameterization. The impact is less during the winter and spring months. The standard deviations for both emissivity representations increase in the spring, possibly due to surface melting, which would impact the horizontal polarization more than the vertical polarization.
4. Selected cases comparing assimilation of a retrieved ice concentration with direct assimilation of brightness temperatures
a. Spatial distribution of the innovations and analysis increments
In this section the spatial distribution of the innovations and ice analysis increments are examined for a particular day taken from assimilation experiments covering the entire ice season. One of the fundamental differences between assimilating an ice concentration computed using the NT2 algorithm and assimilating brightness temperatures directly is that in the former method the observation values (here retrieved ice concentration) are limited to a specific range, here 0%–100%. Close to the limits of 0% and 100% the variability in the retrieved ice concentration is reduced artificially because of the truncation of the distribution. A recent study by Andersen et al. (2007) in an area of the Arctic where ice concentrations were close to 100% found that when the upper limit of 100% was removed, the mean ice concentration from NT2 was 105% with a standard deviation of approximately 5%. This is compared to a mean concentration of 98.9% and a standard deviation of 1.8% when the ice concentration was restricted to be below 100%. This type of bound does not exist for the assimilation of brightness temperatures, which reflect more accurately the true variability of ice conditions, and generate increments over the entire domain. This can be seen in Fig. 7, which shows predicted minus observed (P − O) ice concentration and predicted minus observed brightness temperatures at 0600 UTC 26 January 2007. It can be seen that, when the ice concentration from NT2 is assimilated, the P − O of ice concentration is high along the ice edge, (Fig. 7a), and the increments are mostly along the ice edge (Fig. 7c). When brightness temperatures are assimilated directly there are large values of P − O both along the ice edge and in the ice pack (Fig. 7b), and there are also more significant ice concentration increments in the ice pack, shown in (Fig. 7d), although it is not known how realistic these increments are.
b. Weather effects
The NT2 algorithm has two different ways of accounting for weather effects. The first method is to apply weather filters to the data. There are two weather filters based on empirically determined threshold values for the gradient ratios GR37V22V and GR37V19V. The weather filters cannot change the ice concentration itself other than setting it to zero, and have the impact of reducing the ice concentration along the ice edge. A second method NT2 uses to account for weather effects that applies over both ice and open water is that it uses a radiative transfer model to calculate lookup tables of the calculated gradient ratio and polarization ratio values based on sample atmospheres, ice concentrations, and ice type (thin, first-year, multiyear, or layered). The algorithm then cycles through the values in the lookup tables and compares the calculated gradient ratio and polarization ratio values with those from the observations. The value of ice concentration for which the difference is the smallest is chosen as the retrieved ice concentration. This method produces improved results over the previous NT algorithm, but on a day-to-day basis toggles unrealistically between different atmospheric states (Markus and Cavalieri 2000).
The NT2 algorithm has been developed and calibrated for the Arctic and Antarctica separately (there are two sets of lookup tables), and may not work as well in more temperate ice zones, or in marginal ice zones (Meier 2005; Markus and Cavalieri 2009). An example showing the difficulty the NT2 algorithm has with severe weather effects is shown in Fig. 8 for the 1800 UTC 12 March 2007 analysis. On this day the cloud liquid water in the background state, shown in Fig. 8a, was very high. This produced spurious ice in the NT2 ice retrieval, which appears here in the analyzed ice concentration for NT2, shown in Fig. 8b. To verify that this is spurious ice, it has been checked that there is no ice in the CIS ice charts at that time and location, and it has also been checked that the ice appears very suddenly, and then disappears. The amount of spurious ice in the analysis when using the RTM, shown in Fig. 8c is much less (around 30%) than that in the NT2 analysis (around 70%–80%).
When the RTM is used it may be possible to reduce the positive ice increment due to the cloud liquid water by relinearizing the observation operator. To investigate if relinearization is needed, recall that the approximation was used. By looking at the difference between H(xa) − H(xb) and the validity of the linear approximation can be evaluated. This difference is shown in Fig. 8d for the horizontal 23-GHz channel. It can be seen that the difference is large in the area where the spurious ice appears, suggesting that relinearizing the observation operator and solving for a new analysis state may be beneficial.
5. Verification of the analyses from 3D-Var experiments over the ice season
In this section results from the three data assimilation experiments discussed in section 2e, which are carried out over the period of 5 December 2006–30 June 2007 are compared. The first experiment assimilates ice concentration calculated using the AMSR-E brightness temperatures and the NT2 algorithm. The second and third experiments assimilate the AMSR-E brightness temperatures directly using the radiative transfer model, with the second experiment using the sea ice emissivity values from Eppler et al. (1992), and the third experiment using sea ice emissivity values from our monthly-varying emissivity parameterization.
a. Total ice concentration
Results from evaluating the total ice concentration analysis in terms of the bias and standard deviation of differences with CIS daily ice charts, calculated using values over the 7-month period are given in Table 3. When brightness temperatures are assimilated directly using the emissivity parameterization the averaged statistics are slightly better than assimilating the ice concentration retrieved using NT2.
The bias and standard deviation of ice concentration differences with the CIS daily charts over the season are shown in Fig. 9a as running 30-day means. The difference between the bias for the experiment that assimilates NT2 and that which assimilates the brightness temperatures directly using the variable emissivity are statistically significant at the 95% level over the entire ice season, while the black lines along the horizontal axis of Fig. 9 denote the days when the difference between the variances of the experiments that assimilate NT2 and the brightness temperatures directly using the variable emissivity are statistically significant at the 95% level. Near the beginning of the experiment the biases seem to indicate higher ice concentrations in the analyses that assimilate the brightness temperatures directly using the emissivity parameterization (dashed line) than those that assimilate the NT2 retrievals (solid line). This is because there are some days with a small amount of weather contamination that contributes to spurious ice of low concentration <5% at this time of year. These points are assigned 0% ice concentration in the NT2 algorithm because of the weather filters. This impacts the verification statistics during December because the verification domain is very small at this time of year, as there is only a little ice along the Labrador coast. The positive bias for the direct assimilation of brightness temperatures with the emissivity parameterization is higher than that without at the beginning of the cycle because the weighting given to the observations for the case with the emissivity parameterization is much higher than that with the constant emissivity during December. In Fig. 6 it can be seen that the standard deviations of the innovations are over twice as large for the constant emissivity case than for the emissivity parameterization. It is these standard deviations that are used in our definition of the observation error covariance matrix. Because there are some weather effects in the observations, and our observations are weighted more heavily when the emissivity parameterization is used, the concentration of spurious ice is slightly higher when using the emissivity parameterization.
Over the winter season (January–March) the assimilation of brightness temperatures using the emissivity parameterization produces the best ice concentration analysis out of the three. During this time ice from the northern part of the domain (which formed earlier in the season) advects southward. Hence, the improvements due to the emissivity parameterization are in part due to the positive impact seen in the predicted minus observed brightness temperature scores earlier in the season (see in particular the month of December in Fig. 6). During April and May the standard deviation of the two experiments that assimilate brightness temperatures become quite similar, probably because the ice thickness parameterization of emissivity is no longer used during these months because of the lack of thin ice. The assimilation of brightness temperatures with constant emissivity retains its negative bias compared to the one with the emissivity parameterization, although both have a greater negative bias than when assimilating NT2. It is interesting to note that while NT2 performs well during this time, this is because the minimization in the NT2 algorithm is choosing the thin ice tie points although, according to the ice–ocean model, there is almost no thin ice in the domain. Choosing the thin ice tie points could lead to a higher ice concentration being estimated by NT2, and a lower bias in the ice concentration difference. In June the assimilation of brightness temperatures using the emissivity parameterization again performs much better than the other two methods. At this time of year the model shows significant meltwater on the ice. The ice retrieval (NT2) is known to generate ice concentration estimated during melt that are biased low because the ice and open water signatures are similar to each other. When brightness temperatures are assimilated with the emissivity parameterization, the emissivity calculated during June is much lower than in the winter (see Fig. 2), reflecting the presence of meltwater on the ice. This leads to an improved estimate of the brightness temperatures and a better ice concentration analysis.
b. Ice thickness
The results from evaluating the ice thickness are only shown for the period up to the end of February as this is when the verification data are the most reliable. This is because the ice charts indicate the stage of development, which can only increase or stay constant. When ice starts to melt, the ice becomes thinner, while the stage of development remains constant. The bias and standard deviations of differences with CIS daily charts calculated over the ice season are given in Table 3. The experiment assimilating brightness temperatures with the constant emissivity has the largest bias and standard deviation of the three, whereas averaged statistics from assimilating NT2 ice concentration and assimilating brightness temperatures using the emissivity parameterization are more similar. The bias and standard deviation of ice thickness differences with the CIS daily charts over the season are shown in Fig. 9b as running 30-day means. Even though the bias and standard deviation of the experiment, which assimilates brightness temperatures with the emissivity constant, are lower over the first part of the ice season, there are fewer data points going into the seasonal statistics at this time of year, which is why the seasonal statistics in Table 3 indicate this is the experiment with the overall highest bias and standard deviation. The difference between the bias and variance for experiments assimilating NT2 and those assimilating the brightness temperatures directly using the variable emissivity were found to be statistically significant at the 95% level over the entire ice season.
c. Verification of atmospheric fields
The direct assimilation of brightness temperatures using a radiative transfer model enables analysis increments of water vapor, wind speed, and cloud liquid water to be calculated. In the present study these increments are not used because the ice–ocean model is not coupled to an atmospheric model. However, it is of interest to know if our analysis is closer to an independent estimate than the background. In this case, the independent atmospheric fields are provided by Remote Sensing Systems using retrievals from AMSR-E brightness temperatures (Wentz and Meissner 2000). The forecasts from GEM, which are available every 3 h, were linearly interpolated to the valid time of the retrieval, which is indicated in the retrieval file. The bias and standard deviation averaged over the 7-month period for the background and the analysis are given in Table 4. It can be seen that the analysis of cloud liquid water and wind speed has a reduced bias compared with using the background state, while the bias for the water vapor is higher. The standard deviation for all variables is very similar whether the comparison is performed with either the background state or analysis state.
In this paper, the direct assimilation of AMSR-E brightness temperatures over an ice-covered ocean is presented and compared with assimilating an ice concentration calculated using a retrieval algorithm (NT2) that uses the AMSR-E brightness temperatures as input. The two methods are found to be very comparable with respect to the seasonally averaged statistics when the assimilation of brightness temperatures is carried out using a monthly varying emissivity parameterization.
The direct assimilation of brightness temperatures has several advantages as compared to using a retrieved ice concentration. By using a radiative transfer model, the current forecast of the atmospheric state can be used directly. This allows for better handling of weather effects as compared to assimilating an ice concentration retrieval that uses tabulated atmospheric effects. Another advantage is that by using a footprint operator to assimilate the observations at their different resolutions, the high- and low-frequency channels provide complementary information at the different resolutions. In contrast, when a retrieved ice concentration is assimilated an additional source of error is introduced due to the mismatch in the footprints of the channels used. The main disadvantage associated with the direct assimilation of brightness temperatures is that the method is very sensitive to the ice emissivity and ice temperature, and our a priori knowledge of these two quantities is poor. The method used here to parameterize the emissivity as a function of ice thickness is very simple, and fairly robust because it does not rely on accurate knowledge of the surface conditions. Future work could include continuing with an empirical parameterization, but using a multivariate regression method. For example, the current method was not able to capture the variability of sea ice emissivity in the early spring (April/May). The radiative transfer model should also be modified to use the ice temperature at the appropriate penetration depth for each frequency, and relinearization of the observation operator should also be investigated. A secondary disadvantage of using the radiative transfer model is that it is more computationally expensive than assimilating an ice concentration retrieval, taking approximately 40–45 iterations to reduce the cost function whereas assimilating the ice concentration retrieval takes approximately 15 iterations. However, in the future, when a coupled atmospheric–ice–ocean model is available, this extra cost may be more easily justified if the increments made to the atmospheric variables in the vicinity of ice have a positive impact on weather forecasts. To take full advantage of this, a more sophisticated RTM, like the Radiative Transfer for the TIROS Operational Vertical Sounder (RRTOV; Saunders et al. 1999), should be used after being adapted to include sea ice.
The method used here to calculate the ice emissivity used an estimate of the ice concentration determined by assimilating NT2 ice concentration retrievals. This is however, not necessary, and the same information from assimilating brightness temperatures over a past period could be used instead, although an initial estimate of the emissivity would still be needed. This would make the method independent of NT2, and would allow the observations being used to estimate the emissivity to be independent of those being assimilated. This has been tested for a 1-month interval where forecasts from assimilating brightness temperatures over the past week were used with observations to calculate the emissivity and innovation statistics. Results were found to be similar to the present approach.
The authors thank Leif Toudal Pedersen at the Danish Meteorological Institute for providing the radiative transfer model, and Gregory Smith and the reviewers for helpful comments regarding an earlier version of the paper. This work has been supported by the Canadian Space Agency Government Related Initiative Program (CSA-GRIP).