Baseline Algorithm Definition#
Introduction#
In the following the main algorithm for the multi parameter retrieval is described in detail. The algorithm is based on the works of [Pedersen, 1991, Scarlat, 2018, Scarlat et al., 2017, Scarlat et al., 2020]. The algorithm is divided into several steps, which are described in the following.
Retrieval method#
The retrieval is following a typical scheme with the objective to minimize
where
i.e., \(\mathbf{y}\) is the vector of input brightness temperatures (measurement vector) and \(\mathbf x\) is the vector with the corresponding physical quantities (state vector), \(\mathbf S_e\) is the error covariance matrix of the input brightness temperatures, \(\mathbf S_a\) is the covariance matrix of the a priori values, \(\textbf x_a\) is the a priori value, and \(F\) is the forward operator, which is the compositional forward model described in Forward Model. The optimal solution \(\hat {\mathbf y}\) which minimizes equation (1) (and maximizes the posterior probability density function) is found by iteration. The uncertainty of \(\mathbf y\) can then be obtained as
with \(\mathbf {\hat M}\) is the Jacobian of the optimal solution of (1). The uncertainty of the individual parameters is then given by the diagonal elements of \(\mathbf{\hat S}_a\).
Note
\(\mathbf S_e\) is often set higher than the pure radiometric uncertainty of the brightness temperatures, to account for the uncertainty of the forward model. This can be named an effective error covariance of the measurements. The effective contribution from the forward model error is not evaluated yet, but estimated to be below 2 K for all channels.
Forward Model#
The forward model is the compositional forward model. It consists of individual components, namely the ocean surface, the sea ice, and the atmosphere. At the low frequency channels of the CIMR satellite, the sensitivity to atmospheric parameters is relatively small, but nevertheless the atmosphere needs to be considered. The forward model for ocean and atmosphere for frequencies 6.9, 10.7, 18.7, and 36.5 GHz is used from the AMSR2 ATBD from [Wentz and Meissner, 2000]. The surface contribution to the brightness temperature is given by
where \(C_{\text{ow}}\), \(C_{\text{fyi}}\), and \(C_{\text{myi}}\) are the area fraction of ocean water, first-year ice, and multiyear ice, respectively. \(ε_{\text{ow}}\), \(ε_{\text{fyi}}\), and \(ε_{\text{myi}}\) are the emissivity of ocean water, first-year ice, and multiyear ice, respectively. \(T_{\text{ow}}\), \(T_{\text{fyi}}\), and \(T_{\text{myi}}\) are the brightness temperature of ocean water, first-year ice, and multiyear ice, respectively. \(C_{\text{ow}}\), \(C_{\text{fyi}}\), and \(C_{\text{myi}}\) are adding up to one.
Sea ice type#
For sea ice, the frequency-dependent emissivities for first-year ice and multiyear ice are derived from [Mathew et al., 2009] as
with \(U_{t,h}\) being the temperature-corrected upwelling brightness temperature for the polarization \(p\) at the air temperature at the surface \(T_{C}\) (in °C), \(a_{t}\) and \(b_{t}\) are the frequency dependent coefficients from [Mathew et al., 2009] (see Table 13), and \(ε_{t,p}\) is the frequency-dependent emissivity for the polarization \(p\) and ice type \(t\) from Table 11. [Mathew et al., 2009] did not provide parametrization for L-band, but to match a temperature dependence of emitted radiation at L-band we derive the effective temperature as with (5) using \(a_t=0.1\) and \(b_t=0\) for both, first-year and multiyear ice. The emissivities used for L-band are taken from [Scarlat et al., 2020] with slight modifications and are given in Table 11.
Note
\(T_{C}\) in this parametrization is the IST in the retrieval as the only temperature dependence. The assumption is that in thermal equilibrium the IST is equal to the air temperature at the surface.
Sea ice thickness#
The ice thickness dependence at all frequencies is introduced via modification of the emission from the ice surface for the first-year ice fraction. The ice emissivity is modified by the ice thickness \(\text{SIT}\) according to
with the index \(p\) indicate polarization (\(h\) or \(v\)) the coefficients \(a_{p}\), \(b_{p}\), and \(c_{p}\) from [Scarlat et al., 2020] (see Table 12). To combine ice temperature and ice thickness dependence, the ice emissivity is modified with the ice thickness by substituting \(a_{p}\) with the ice temperature sensitivity \(U_{T, t, p}\) from (5) for \(t=\text{FYI}\). This was not performed in [Scarlat et al., 2020] but is essential for the minimization of the cost function to not introduce discontinuities in the forward model. The MYI emissivity is not affected by the ice thickness in this forward model, as it was not part of the thickness sensitivity study in [Scarlat et al., 2020]. As a consequence, the retrieval of ice thickness accounts only for the FYI fraction and thus is very noisy when the FYI fraction is small. However, the winter season during freeze-up conditions, this is not much of a restrictions as most of the ice is first-year ice thin multiyear ice is scarce.
Ocean surface model#
While for the sea ice an empirical model is used for the emissivity, for the emission from the ocean the model uses the Fresnel reflection coefficient as a basis, which relies on the dielectric constant of the sea water. The emission from calm sea water after [Meissner and Wentz, 2012] is given by
with \(ε\) being the dielectric constant of the sea water. To account for the roughness and other disturbances on the ocean surface, the power reflectivity at each polarization \(R_{0p}=|r_p|^2\) needs to be adjusted. The contribution of different ocean surface types are modeled by different parameters again following [Wentz and Meissner, 2000]. The reflectivity is a composition of foam covered ocean and clear ocean. With a reduction of the ocean surface emission through the foam by a factor κ, the composite reflectivity is given by
with \(f_{\text{foam}}\) being the fraction of the ocean surface covered by foam and \(R_{\text{clear}}\) being the reflectivity of the clear ocean. With a small loss from diffraction expressed in the term \(β\) we can express the reflectivity as
with \(R_{\text{geo}}\) being the reflectivity from a standard geometric optics model[Wentz, 1975] defined later in (10). The combination can then be expressed by combining the foam and the diffraction term into one quantity \(F=f_{\text{foam}} + \beta -f_{foam}\cdot β - f_{foam} κβ\), which is a monotonic function of wind speed and is addressed by [Wentz and Meissner, 2000] as catch-all term. They determined \(F\) empirically from experiments with various radiometers. A fit for F is given by
a quadratic spline with knots at \(W_1=3~\text{m/s}\) and \(W_2=12~\text{m/s}\) for v- and \(W_1=7~\text{m/s}\) and \(W_2=12~\text{m/s}\) for h-polarization. The coefficients \(m_1\) and \(m_2\) are given in Table 10 in the appendix.
The geometrical optics surface roughness (\(R_{geo}\) from (8) is used as modeled by [Wentz and Meissner, 2000] as
Here \(R_0\) is the specular reflection, θ\(_i\) is the incidence angle, T\(_S\) is the sea surface temperature, and W is the wind speed. The coefficients for each polarization and frequency are in the appendix in Table 9.
Dielectric constant of sea water#
The dielectric constant of sea water depends on salinity and temperature [Meissner and Wentz, 2004] as
with ε\(_1\) as the intermediate frequency dielectric constant, ε\(_S\) as the static dielectric constant, ε\(_\infty\) as the high-frequency dielectric constant, \(\nu_1\) and \(\nu_2\) as first and second Debye relaxation frequencies, respectively and σ as the conductivity. Fits for S=0 are given by [Meissner and Wentz, 2004] as
where the salinity dependence is modeled [Meissner and Wentz, 2004] as
However, [Meissner and Wentz, 2012] published an update to the salinity dependence of ε\(_S\) and ν\(_1\) as
We use this most recent update in the forward model.
For conductivity we follow [Meissner and Wentz, 2004] and [Meissner and Wentz, 2012] who used the earlier work by [Stogryn et al., 1995] as
with
Atmosphere#
To get back to emissivity in order to account for the atmospheric contribution to the brightness temperatures at surface level, the brightness temperature is just devided by the ice and ocean surface temperatures
With \(C_{\text{myi}}+C_{\text{fyi}}=\text{SIC}\) and \(\text{SIC}\cdot C_\text{myi} = \text{MYI}\) being part of the state vector (2), the state defines the surface area fraction of all three considered surface types.
The reflectivity of the surface, \(R_{\text{surf}}\), is given by
The atmospheric contribution to the brightness temperature is calculated from the parametrization from [Wentz and Meissner, 2000]. They fitted the downwelling and upwelling effective temperature using a least-squares fit to the TWV. The fit is given by
where \(T_v = 273.16+0.8337 V - 3.029\cdot 10^{-5}V^{3.33}\) for V<48 and \(T_v=301.16\) for V>48, \(\zeta(x)=1.05x(1-x^2)/1200\) for \(|x|<20\ \text{K}\) and \(\zeta(x)=\text{sign}(x)\cdot 14\ \text{K}\) for \(|x|>20\ \text{K}\).
The absorption by oxygen is given by
The vapor absorption is given by
The liquid water absorption is given by
with \(L\) being the liquid water path.
The total atmospheric attenuation is given by a combination of the individual terms from (20), (21), and (22):
with \(\theta\) being the incidence angle.
An appropriate approximation of the atmospheric integration for up- and downwelling brightness temperature used by [Wentz and Meissner, 2000] is given by
This makes the upwelling and downwelling brightness temperature differ only slightly by some Kelvin.
Addition of the L-band forward model#
The addition for 1.4 GHz was done by [Scarlat et al., 2020] and is based on [Ruf, 2003].
The atmospheric attenuation for L-band is
with \(V\) being the TWV in mm and \(\theta\) being the incidence angle.
The up- and down-welling brightness temperature for L-band is given by
with ST being the surface temperature in K. Over open ocean, this is SST and over ice it is IST. The composite upwelling and down-welling brightness temperature is then calculated with the weight of the ice concentration.
The effect of wind roughening for L-band over open ocean is given by
with \(u\) being the wind speed in m/s, \(\theta\) the incidence angle in degrees, and \(E_{0,h}\) and \(E_{0,v}\) being the emissivity of the surface in the horizontal and vertical polarization from Eq. (7).
Brightness temperature at instrument level#
Bringing together the different models for all frequencies, the brightness temperature at the instrument level is given by
with \(p\) being the polarization and \(T_c\) being the cosmic background temperature. \(T_{b,s}\) is the surface emitted brightness temperature for the combined surface and \(\epsilon_p\) is the emissivity of the combined surface for the polarization \(p\).
Correction for variation in incidence angle#
While the algorithm is designed for a fixed incidence angle, the incidence angle is varying slightly depending on the feed horn which is used for the observation which, in turn, changes the brightness temperatures. We attempt a straightforward correction for the incidence angle dependence of the brightness temperatures by estimating the surface reflectivity using the polarization information for a given observation. The assumption here is a simple Fresnel reflection model without an atmospheric contribution. First, we define the Fresnel reflection coefficient of the surface with given effective permittivity \(\varepsilon\) and incidence angle \(\theta_i\) as described in (7). The effective permittivity is used here which is a simplification calculated with the assumption that the atmosphere is transparent. We use an effective surface temperature \(T_{\text{eff}}\) which we set to \(T_{\text{eff}}=240\) if \(T_{b,v}<220\) K otherwise it is set to \(T_{\text{eff}}=T_{b,v}+20\) K as an approximation. Then we have
which we can solve for a given \(θ_i\) to get the effective permittivity as
The corrected brightness temperature for the incidence angle \(θ_{\text{ref}}\) at vertical polarization is then given by
with \(θ\) being the incidence angle of the observation. In theory, we expect the correction factor \(ζ\) to be equal to the effective temperature \(T_{\text{eff}}\). Our only source for the incidence angle dependence is a simulated scene from the SCEPS project, referred to as the SCEPS polar scene in this document. In this scene, the small differences in incidence angle between different scans are included and modeled in the brightness temperatures. For \(T_{b,v}\) this correction works for all CIMR frequencies on the SCEPS polar scene with a correction factor \(ζ=T_{\text{eff}}\) as expected. For \(T_{b,h}\) a slightly different correction factor was found in order to remove the incidence angle dependence of the brightness temperatures in the SCEPS polar scene. The origin of this factor is not yet understood and is subject to further investigation. In this version of the algorithm, the correction factor \(ζ\) for \(T_{b,h}\) is shown in the table Table 1 and is applied to the \(T_{b,h}\) for each channel.
Band |
Correction factor \(ζ\) |
---|---|
L_BAND |
\(T_{\text{eff}}\) |
C_BAND |
\(T_{\text{eff}} - T_{b,h}(θ)\) |
X_BAND |
\((T_{\text{eff}} - T_{b,h}(θ))/0.9\) |
KU_BAND |
\((T_{\text{eff}} - T_{b,h}(θ))/1.4\) |
KA_BAND |
\((T_{\text{eff}} - T_{b,h}(θ))/1.8\) |
This method requires a non-linear optimization for each observation, which can be computationally costly. It is expected to work without modification only on the SCEPS simulated scene as it is probably similar to the method used to create the incidence angle dependence in the scene. A real CIMR scene might require a different approach.
Note
The proposed method here is solely to correct the incidence angle dependence of the brightness temperatures, the involved effective permittivity \(ε_{\text{eff}}\) and effective temperature \(T_{\text{eff}}\) should not be confused with the physical permittivity or the temperature of the surface, respectively.
CIMR Level-1b re-sampling approach#
The CIMR Level-1b data is first resampled to the C-band footprints using nearest neighbor resampling. The target footprint with C-band is a compromise between resolution and accuracy given sensitivities of the inversion of the forward model. The main restraint being the L-band channel which has a high influence on the retrieval, on ocean because of SSS and sea ice due to SIT. Higher up sampling of L-band than to the C-band resolution would require accurate knowledge about the antenna pattern of the instrument and a trade-off between accuracy and resolution.
Level-2 end to end algorithm functional flow diagram#
The diagram below shows the functional flow of the algorithm.
Functional description of each algorithm step#
Input Data: The algorithm begins by gathering input data from the CIMR Level-1b products, which include top-of-atmosphere (TOA) brightness temperatures (TBs) and their associated errors (TBe). Additionally, external data sources such as ECMWF analysis and historical ERA5 datasets are incorporated to provide necessary atmospheric and environmental context.
Resampling: The input data is resampled to the C-band footprints using a nearest neighbor approach. This step ensures that the data is aligned with the target resolution, balancing the need for accuracy and resolution, particularly for the L-band channel which significantly influences the retrieval process.
Forward Model: The forward model simulates the expected brightness temperatures based on the input data and the current state vector. It incorporates atmospheric attenuation, surface emissivity, and other relevant physical processes to generate a consistent state of the observed brightness temperatures.
Cost Function: The cost function quantifies the difference between the observed brightness temperatures and those predicted by the forward model. It incorporates measurement uncertainties and a priori information to evaluate the quality of the fit. This step is crucial for optimizing the state vector to minimize discrepancies between observed and modeled data.
Optimal State: The optimization process iteratively adjusts the state vector to minimize the cost function. This involves using the Levenberg–Marquardt algorithm to find the best-fitting parameters that explain the observed brightness temperatures. The resulting optimal state vector represents the most accurate estimate of the geophysical variables.
Output Data Level 2: Once the optimal state vector is determined, the algorithm generates output data, which includes geophysical variables derived from the brightness temperatures, uncertainties associated with these variables, residuals of the brightness temperatures, and quality flags indicating the reliability of the retrievals.
Algorithm assumptions and simplifications#
Compared to individual state-of-the-art algorithms for the retrieval of individual quantities, the multi parameter retrieval uses a simplified approach. Particular trade-offs include:
Some ocean parameters contain many empirical values, developed and validated for different instruments which might have to be adjusted to match CIMR instrument characteristics.
The sea ice thickness is retrieved as a combination of first-year ice emissivity and the ocean emissivity. As this, it is prone to noise in the ocean emission, which can lead to erroneous occurrence of sea ice of low thickness, as the first-year ice emissivity is identical to the ocean emissivity at 0 cm thickness.
The empirical parametrization of sea ice thickness stems from uncertain atmospheric conditions, so that the ice thickness dependence of higher frequency channels might be confounded with the dependence on the atmospheric conditions.
The sensitivity of the CIMR channels to CLW, in particular over sea ice is low, so that the uncertainty of CLW is high. In addition, the effect of decreased CLW over first-year ice is similar to the effect of an increased MYI fraction, so that these two signals cannot be adequately separated in the current retrieval.
In the present version of the algorithm, the forward model is assumed near-linear in its characteristic and the retrieval is unconstrained. This might lead to a bias in the retrieval of the parameters when other parameters might converge to nonphysical values. In test cases, this happened mostly with
IST > 273.15 K
SST < 273.15 K
TWV < 0 kg/m\(^2\)
There are more modern approaches for the forward model of the ocean emission, in particular at L-band. This will influence the retrieval of the ocean parameters, in particular SSS and SST. The current algorithm is the approach of [Ruf, 2003] and [Scarlat et al., 2020], but a switch to another parametrization like [Meissner et al., 2018] can be considered. However, the SSS is the parameter with the weakest signal in the forward model, so that other parametrization have to be improved, to justify this effort.
The forward model is restricted to winter conditions over sea ice. The summer conditions are highly variable across frequencies and cannot be adequately described by the current forward model. In particular the definition of first-year ice and multiyear ice is ambiguous in melting conditions. In addition, the SIT parametrization is only valid for winter conditions.