Modeling Snow and Ice Microwave Emissions in the Arctic for a Multi‐Parameter Retrieval of Surface and Atmospheric Variables From Microwave Radiometer Satellite Data

Monitoring surface and atmospheric parameters—like water vapor—is challenging in the Arctic, despite the daily Arctic‐wide coverage of spaceborne microwave radiometer data. This is mainly due to the difficulties in characterizing the sea ice surface emission: sea ice and snow microwave emission is high and highly variable. There are very few data sets combining relevant in situ measurements with co‐located remote sensing data, which further complicates the development of accurate retrieval algorithms. Here, we present a multi‐parameter retrieval based on the inversion of a forward model for both, atmosphere and surface, for non‐melting conditions. The model consists of a layered microwave emission model of snow and ice. Since snow scattering and emission effects, as well as temperature gradients, are taken into account, a high variability in brightness temperatures can be simulated. For ocean regions and the atmosphere existing parameterized forward models are used. By using optimal estimation, the forward model can be inverted allowing for the simultaneous and consistent retrieval of nine variables: integrated water vapor, liquid water path, sea ice concentration, multi‐year ice fraction, snow depth, snow‐ice interface temperature and snow‐air interface temperature as well as sea‐surface temperature and wind speed (over open ocean). In addition, the method provides retrieval uncertainty estimates for each retrieved parameter. To evaluate the forward model as well as the retrieval, we use the extensive data sets acquired during the year‐long Arctic expedition Multidisciplinary drifting Observatory for the Study of Arctic Climate (2019–2020) as a reference.

and lapse rate feedback) but is hampered by the sparseness of in situ measurements and the challenges using satellite microwave observations for atmospheric sounding over sea ice (Crewell et al., 2021;Kang et al., 2023).The surface contributions to the signal dominate the microwave brightness temperatures measured at the top of the atmosphere compared to the atmospheric emissions, with the surface signal being transmitted through the atmosphere even up to 480 GHz at the poles because of the dry atmosphere (Wang et al., 2017).Microwave-based retrievals of geophysical parameters in the Arctic are challenged by the variability of emissions of both, snow and sea ice.Often, empirical parameterizations or climatologies for the surface emissions are used (e.g., Kilic et al., 2020;Wang et al., 2017).While the strong surface signal is demanding for atmospheric parameter retrievals like for water vapor, it can prove useful regarding information about variables like snow depth (SND) on sea ice, that play important roles in feedback mechanisms, e.g., the insulating role of snow affecting sea ice growth.Thus, an integrated approach, retrieving both atmospheric and surface variables, is taken in this study.
We build upon works of Scarlat et al. (2017Scarlat et al. ( , 2020) ) who took results from Mathew et al. (2009) to derive an empirical and spatially non varying, effective emissivity and emitting layer temperature of the sea ice for the surface component of a forward model.This model is inverted by an optimal estimation method to simultaneously retrieve several geophysical parameters including sea ice concentration (SIC), that is, the percentage of an ocean area covered by sea ice, and total water vapor (TWV), that is, the total amount of water vapor in an atmospheric column, from brightness temperatures measured at six frequencies ranging from 6.9 to 89 GHz.When studying the derived atmospheric quantities from this method, we note a bias for water vapor when compared to radiosonde measurements during field campaigns over sea ice, which is strongly reduced over open ocean (Crewell et al., 2021).Also, this approach shows unrealistic high values of liquid water path (LWP) over first-year ice (FYI), and thus a false gradient when the ice type changes from predominately first-year to multiyear ice as we will see later in this manuscript.This is due to the assumed fixed surface emissivities (that are only distinguished by ice type), especially at the high-frequency channels, which carry most of the information about the atmospheric state.At these frequencies in particular, both scattering and emission in the snow play a role, but also the vertical temperature profiles in both snow and ice that depend on the depth of the insulating snow.In order to better represent these highly variable surface emissions in the forward model, we exchange the empirical parameterization with the Microwave Emission Model of Layered Snowpacks extended to sea ice, called MEMLS_ice model (Mätzler & Wiesmann, 1999;Tonboe, 2010;Tonboe et al., 2006;Wiesmann & Mätzler, 1999) and include SND and snow-ice interface temperature as new retrieval parameters.The MEMLS_ice emission model uses a stack of planar homogeneous layers characterized by the physical snow and sea ice quantities relevant for scattering and absorption, like the layer thickness (m), correlation length (mm), a measure of scatterer size and distribution, density (kg m −3 ), salinity (ppt), and temperature (K).
Here, we present first results of this physical approach for integrated retrievals of surface and atmosphere parameters for the cold season over the Arctic Ocean (October to April).First, an idealized layer input to MEMLS_ice is established based on literature and campaign data for the layer parameters.The sensitivity to the model layer choices is analyzed and used as uncertainty estimate in the inversion method.To showcase the impact of the physical ice emission model we make use of the extensive data sets collected during the year-long Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition (Nicolaus et al., 2022).Here, we study the brightness temperatures simulated by the new forward model using ground truth input from MOSAiC and compare them to satellite observations.Then, the model is inverted and the retrieval output is analyzed.Finally, we apply the retrieval on Arctic-wide available satellite data to generate daily maps of all retrieved quantities.

AMSR2 Satellite Data
We use brightness temperatures measured by AMSR2, the microwave scanning radiometer on the GCOM-W1 spacecraft from the Japan Aerospace Exploration Agency (JAXA), launched 18 May 2012.The dual-polarized microwave radiometer has an instantaneous fields of view (IFOV) ranging from 62 × 35 km at 6.9 GHz to 5 × 3 km at 89 GHz and a swath width of 1,450 km.Its near-polar, sun-synchronous sub-recurrent circular orbit is at an altitude of about 700 km.
We use the swath data, both ascending and descending, of the Level 1R product matched to the resolution of 6.9 GHz (Maeda et al., 2016), corresponding to an IFOV of 62 × 35 km.This product is resampled using the 10.1029/2023EA003177 3 of 29 antenna gain patterns for each channel so that the brightness temperatures for all frequencies have the same field of view.We use the land mask contained in the product to mask out land observations.In this analysis, we also use data of overpasses over the icebreaker Polarstern (Knust, 2017) during MOSAiC.For every overpass, the satellite measurement closest to the vessel's hourly position is taken and used in Section 4.1.There are up to seven overpasses of Polarstern per day.Additionally, we use daily gridded Level 1B data projected onto a polar stereographic grid with a nominal resolution of 25 km to derive a climatology in Section 3.2.2.

ERA5 Reanalysis Data
From the European Centre for Medium-Range Weather Forecasts fifth generation reanalysis ERA5 (Hersbach et al., 2020) we obtain climatological mean values further described in Section 3.2.2.

Merged Warren-AMSR2 Climatology: Snow Depth
For a climatology of SND, we use a merged product described in Hendricks and Paul (2022) of monthly SND data based on merging the Warren snow climatology (Warren et al., 1999) and daily SND over first-year sea ice from AMSR2 data (Rostosky et al., 2018).

MOSAiC Campaign Data
The MOSAiC expedition (Nicolaus et al., 2022;Rabe et al., 2022;Shupe et al., 2022) provided the unique opportunity of measurements in the central Arctic in winter: the icebreaker Polarstern drifted passively with the packice anchored to a second-year ice floe for a full year from October 2019 to September 2020, with one relocation during the summer.The main site for data collection was the Central Observatory (CO) within a few kilometers of the ship.Predominantly, this was second-year ice (Krumpen et al., 2021).The data sets used in this study as "ground truth" are described in the following.Note that they have different spatial extents, but we assume that the data is representative of the satellite footprint.Therefore, we also use estimations of uncertainties in Section 4.1.

Humidity and Temperature Profiler (HATPRO) Radiometer: TWV and LWP
We use TWV and LWP retrieved from the ground-based low frequency HATPRO microwave radiometer operated onboard the research vessel Polarstern during MOSAiC.More information on this data can be found in Walbröl et al. (2022).Here, we are using only data that was not flagged as bad data in the quality check.

Snow and Ice Mass Balance Apparatus (SIMBA) Buoys: SND and Snow-Ice Interface Temperature
For SND and snow-ice interface temperature measurements during the MOSAiC campaign we use data from 19 SIMBA buoys, which are thermistor string type ice mass balance buoys that were deployed in proximity to Polarstern (Distributed Network [DN] and the CO) (Lei et al., 2021a).Note that not all buoys were operational at the same time but that we have data from at least 10 buoys at different locations within the CO at the same time during our period of interest.

TerraSAR-X ScanSAR Images: SIC and Multiyear Ice Fraction
We use 49 classified scenes for November 2019-March 2020 from Guo et al. (2023) based on X-band TerraSAR-X scanning synthetic aperture radar images in a 70 × 70 km square around the CO.We resampled the data to daily values by linear interpolation and fill missing values with the mean of the whole time series.The SIC is estimated using the fraction of classified lead ice (representing openings, i.e., an upper bound for the open water fraction), while the multiyear ice fraction is estimated to be the sum of the fractions of the deformed and heavily deformed ice classes.Measurements of WSP and SST are provided by the vessel's synoptical reports (Schmithüsen, Raeke, & Kieser, 2021;Schmithüsen, Rohleder, & Hausen, 2021;Schmithüsen, Schröter, & Wenzel, 2021).

Infrared Thermometer: Surface Temperature
For the snow-air interface temperature we use skin temperature calculated based on infrared thermometer data measured at the Met City location within the CO and from three additional stations in the DN (with average distances to the vessel between 10 and 23 km) (C.Cox et al., 2023aCox et al., , 2023bCox et al., , 2023cCox et al., , 2023d;;Herrmannsdörfer et al., 2023).

N-ICE 2015 Campaign Data
In addition to the MOSAiC measurements we use radiosonde data from the Norwegian young sea ICE cruise (N-ICE 2015) that was conducted north of Svalbard in the Atlantic sector of the Arctic in thin, first-year sea ice from January to June 2015 (Cohen et al., 2017;Hudson et al., 2017).The TWV is calculated from the radiosonde data using the saturation vapor pressure formulation by Hyland and Wexler (1983) and we restrict our analysis to January until April.

Round Robin Data Package
The Round Robin Data Package (RRDP, Pedersen et al., 2019) is a collection of co-located sea ice measurements from buoys, flight campaigns, melt pond satellite retrievals, 100% and 0% SIC reference areas, satellite microwave brightness temperatures and backscatter, as well as atmospheric reanalysis data (ERA5 and its predecessor ERA Interim), produced and maintained by ESA's Climate Change Initiative project.It also includes snow and ice data from NASA's Operation Ice Bridge campaigns (Kurtz et al., 2013).We use the RRDP Version 3 and constrain our analysis to the months October-April.

Methodology
For retrieving physical state parameters expressed as x from satellite-measured brightness temperatures TB, we invert the forward model (1) using an optimal estimation method.In this section, forward model and inversion scheme are described.Model and inversion are implemented using the programming language Julia (Bezanson et al., 2017) allowing for fast computations needed in satellite retrievals.

Forward Model
A satellite observes upwelling brightness temperatures, TB, at an incidence angle θ for a given frequency f at the top of the atmosphere given by (Wentz & Meissner, 2000) where TB U is the upwelling emission of the atmosphere, τ is the total transmittance, TB Ω is the downwelling radiation scattered upward by the surface, and TB surf. is the radiation emitted by the surface.These expressions allow us to separate the terms into an atmospheric part (TB U and τ), a surface contribution TB surf.and a term containing both as TB Ω contains atmospheric downwelling radiation as well as the surface reflectivity.
We consider the surface contribution as a linear combination of sea ice and open ocean (OW) weighted by SIC as also done in Scarlat et al. (2017Scarlat et al. ( , 2020)): with TB ice being the brightness temperatures of the sea ice surface and TB OW being the brightness temperatures of the ocean surface.Both terms include the reflected downwelling radiation of the atmosphere by the respective surfaces.We model the atmosphere and ocean following Wentz and Meissner (2000) with recent improvements (Meissner & Wentz, 2004, 2012).The surface temperature is given by the weighted mixture of SST and snow-air interface temperature.We now focus on modeling TB ice in the following Section.

Surface Model
Our aim is to construct a surface model setup that (a) covers the most important quantities influencing the brightness temperatures in order to simulate plausible brightness temperatures, polarization differences and their variabilities and (b) reflects a realistic snowpack including its temperature gradients by using model parameters from literature.Our chosen setup is characterized by four input variables: the snow-air interface temperature T sa , the snow-ice interface temperature T si , the ice type (either FYI or multiyear ice [MYI]) and the SND.We use the MEMLS_ice model to simulate brightness temperatures at the top of multiple stacked planar layers of ice and snow.MEMLS_ice calculates the radiative transfer in a multi-layer media, including internal reflections (Fresnel equations) and scattering.Absorption and scattering coefficients are hereby derived from six-flux theory, taking all space directions into account.The model is described in Wiesmann and Mätzler (1999) and Mätzler and Wiesmann (1999) and was extended to include sea ice by Tonboe et al. (2006) and Tonboe (2010).
The model input consists of downwelling brightness temperatures that are obtained from the atmospheric model (see above) and parameters for each layer, namely thickness, density, correlation length, liquid water content, salinity, and temperature.The incidence angle is set to the satellite incidence angle of 55°.Snow and ice are highly variable both temporally, for example, due to snow metamorphism or ice growth, and spatially, for example, because of different ice types or horizontal in-homogeneity of the snow surface.In our retrieval method we have not enough degrees of freedom in our satellite measurements to retrieve all snow and ice parameters that are needed as input for the MEMLS_ice model for it to produce a realistic output.Therefore we constrain some of them.Using relations between these parameters (i.e., their co-variances) to reduce their effective number is hindered by the high variability of the atmosphere that partly determines the ice and snow structure during initial formation (Fuhrhop et al., 1998).It is therefore clear that our layer setup represents a simplified, idealized ice and snowpack.
We model snow and ice using four and five layers, respectively, in order to capture temperature and salinity profiles as well as layering observed in snow on Arctic sea ice (King et al., 2020;Merkouriadi et al., 2017).This is especially important as the radiation from the surface is the result of emission, scattering and reflection within the snow and ice pack, so that each layer contributes to the observed brightness temperatures (Hallikainen & Winebrenner, 1992).For different frequencies the relative contributions vary and can be related to "penetration depths": low frequencies have larger penetration depths than high frequencies influencing TB differences among them.
All fixed model parameters are shown in Table 1 and are chosen from literature or observations as discussed below.To estimate the influence of these model parameters on the model output (the top of atmosphere brightness temperatures) and to quantify the uncertainty that we are introducing by not being able to account for the model parameters' natural variability, we perform model simulations which are described in Section 3.2.5.There we measure the TB sensitivity to the model parameters, which can be later used in the optimal estimation scheme.

Temperatures
The snow air interface temperature and the snow-ice interface temperature are free parameters and are retrieved in the final satellite product.They are treated as independent parameters, that is, we are not assuming thermal equilibrium within the snow and ice.This allows us to model changing satellite brightness temperatures due to temperature changes in the upper snow layers that do not affect the ice layer's temperatures, which is of particular importance for the higher frequency channel TBs.The temperature profile in the snow is assumed to be linear.Thus temperature changes leading to non-linear temperature profiles cannot be properly modeled with our approach.The same is true for the ice temperature profile, where we assume a temperature of 271.35K for the ocean, the freezing point of seawater with a salinity around 35 ppt.

Snow Parameters
While SND is a free parameter, the other parameters related to the microphysical structure of snow, namely exponential correlation length and density, are chosen in a way to resemble a layer of wind slab snow on top, a faceted layer in the middle and a layer of depth hoar snow at the snow-ice interface, as it was observed on Arctic sea ice (King et al., 2020;Merkouriadi et al., 2017).Note that in reality the distributions and relative thicknesses might vary depending on the ice type and so does the density.For example, higher bulk density of snow over second-year ice compared to FYI has been observed (Merkouriadi et al., 2017).For simplicity and to avoid over-emphasizing the ice type parameter, such a distinction is not made and the snowpack model parameters are independent of the ice type.For snow densities, we take average values from field measurements in April (King et al., 2020).We note the high variability in both the density distributions and fractions to the overall SND (e.g., King et al. (2020) find that the depth hoar layer for MYI makes up 50% of the snowpack but only 20% over FYI).We simplify by setting the fractions of the three layer to one third, and we take average values of FYI and MYI.In our model, wind slab and faceted snow has the highest density (327 and 314 kg m −3 ) and depth hoar the lowest (250 kg m −3 ), yielding a bulk density of 297 kg m −3 in the range of literature values (e.g., Warren et al., 1999).The snow permittivity of each layer in the model, as a function of snow density, is calculated as a dielectric mixture of air and ice described in Wiesmann and Mätzler (1999) following Mätzler (1998).Over FYI, brine-wetted snow has been observed (Barber et al., 1995;Crocker, 1992;Geldsetzer et al., 2009) and its effect on satellite altimetry was discussed in Nandan et al. (2017).The highest salinities are found in the bottom 4-8 cm of the snow cover and range from 1 to over 20 ppt (Barber et al., 1995;Geldsetzer et al., 2009).As our implementation of the MEMLS_ice model does not include saline snow permittivity models, we implemented a dielectric mixture model assuming ellipsoid brine inclusions in the snow.We choose an upper limit of 6 cm for the saline snow layer over FYI with a salinity of 11.4 ppt.A discussion of this can be found in Appendix A. Over multiyear ice, we assume no brine in the snow.
The scattering in the snow is strongly depended on the microstructure of the snow and MEMLS_ice uses the exponential correlation length as corresponding parameter.Within MEMLS_ice, different scattering mechanisms are parameterized.We model the scatterers of the upper two layers as spheres and for the bottom layer, which we want to resemble depth hoar, we choose shell scattering, described in Wiesmann and Mätzler (1999).
The correlation lengths are estimated based on SnowMicroPen MOSAiC measurements from October to April (Macfarlane et al., 2021).Snow surface and snow bottom at the ice interface were automatically detected in the SnowMicroPen measurements and then the snowpack was simply divided in three layers of equal thickness.Mean correlation lengths for each layer were computed using the parameterization given in Proksch et al. (2015) and were then averaged over the whole time series.The values fit to the ones given by Mätzler (2002) for Alpine snow for the faceted and depth hoar snow type.

Ice Parameters
Solid ice, liquid brine, and air bubbles make up the sea ice.As that composition changes over time, for example, due to desalination as result of expulsion of brine, first-year and multiyear ice differ in their microwave emission signature.Thus, ice type is the free parameter determining the microphysical parameters that are assumed for the model ice layers.We are using two ice classes: FYI and multiyear ice.
Sea ice thickness defined as the distance from ice surface to ice bottom, is strongly variable and ranges from a few cm (new ice/nilas) over 1-3 m for level first-and multiyear ice, to tens of meters for ridged ice.Here, we assume a 1 m ice layer for FYI and 2.5 m for multiyear ice.Ridges can be up to tens of meters and individual adjacent floes can strongly vary in their mean thicknesses (Tucker et al., 1992) but we argue that on a satellite footprint it is eligible to assume some mean thickness.The main influence of ice thickness on the modeled TB is given by the thickness-dependent temperature gradients.Due to the inclusion of the snow-ice interface temperature T si as free parameter, this influence is partly reduced.The natural seasonal dependence of sea ice thickness is currently not implemented and the used thicknesses can be understood as some plausible mean winter value for FYI and multiyear ice comparable to Kwok et al. (2020).
For the FYI density, we follow the literature value of around 916 kg m −3 (Alexandrov et al., 2010).Multiyear ice that survived one melting season usually contains more air inclusions and has a lower density above sea level.
The first 25 cm of the ice in our setup can be considered freeboard with lower densities (here the upper ice layer density is constant at 895 kg m −3 ), while the other layers are below sea level with values comparable to FYI (Jutila et al., 2022).Even lower densities for the multiyear ice part above the waterline have been reported (720 to 910 kg m −3 (Timco & Frederking, 1996)).The very low densities values mainly occur in ridges in deformed ice areas, which do not dominate the areal average we are interested in.Thus, we assume that our multiyear ice density value is a good approximation for the satellite footprint scale average.
Saline ice mainly impacts the dielectric loss term.Nakawo and Sinha (1981) observed a quasi-stable value within a few weeks after ice formation with top layer ice salinities of around 11.4 ppt, and we distribute the salinity in the ice with a C-shaped salinity profile model for FYI (Weeks & Lee, 1962).Salinity in multiyear sea ice is much lower than in FYI especially near the surface because of the summer melt processes (G.F. N. Cox & Weeks, 1974).We choose a linear model to distribute salinity in multiyear ice with 1 ppt in the top layer and 4 ppt in the bottom ice layer.
The permittivity of the ice layers are given by the dielectric mixture of ice, air and brine.The permittivity of ice is calculated following Matzler and Wegmuller (1987) with the imaginary part given by Hufford (1991).
The permittivity of brine is calculated following Stogryn and Desargant (1985).The brine volume is calculated as expressed in Stogryn (1987).For FYI, the permittivity of ice and brine is calculated with a dielectric mixture model assuming prolate spheroid inclusions with an axis fraction (i.e., the ratio of the semi-major to the semi-minor axis) of five to resemble brine pockets (Light et al. (2003) also treat brine pockets as prolate ellipsoids).For multiyear ice the effective permittivity is calculated based on effective-medium theory of Polder and van Santeen (1946) assuming spherical inclusions of air in ice.
The scattering mechanisms in ice are estimated by the improved Born approximation and are different depending on the ice type: for FYI small brine pockets are the dominant scatterers, while for multiyear ice air bubbles scatter the radiation following (Tonboe, 2010).Note that the parameterization for ice assumes scattering at spheres and calculates an effective permittivity (for FYI done assuming randomly oriented needles), even though the dielectric mixture models assumes randomly oriented ellipsoids for FYI.The assumed correlation lengths are taken from Rostosky et al. (2020) based on ice core data, where the FYI values are in agreement with values reported by Tucker et al. (1992).

Additional Model Parameters
As our model is designed for the winter season, wetness (caused by meltwater in the ice or snowpack) is set to zero for all layers.Upper ocean salinity is of importance only for frequencies lower than the ones used in this study (Kilic et al., 2021).We use a value of 34 ppt and calculate the ocean permittivity and conductivity following (Meissner & Wentz, 2004).We use the model without taking into account coherence effects.This is necessary to avoid oscillations due to coherent reflections at the parallel layers as soon as one layer gets smaller than 2 cm.Additionally, the roughness of layer interfaces is not considered in our model.In general, we expect that the high frequencies (in this case at 89 GHz) are most challenging to model in the MEMLS_ice model.This is due to the small wavelength which can be of the order of, for example, the brine inclusions or the interface roughness, so that other scattering mechanisms (e.g., geometrical optics) might come into play which are not yet parameterized.

Inversion Scheme: Optimal Estimation Method
Inverting the forward model leads to retrieved parameters that are physically consistent in model space, for example, ice type fractions add up to the total SIC.Here, we use an optimal estimation method for the inversion that additionally retrieves an estimate of the uncertainty for each parameter.The inversion method is a Bayesian approach based on Rodgers (2000) using the conditional probabilities where P(x|y) is the probability of a state x (state vector) given measurements y (measurement vector, here: brightness temperatures of different frequencies and polarizations), P(y|x) is the conditional probability density of y given x, and P(y) and P(x) are the prior probability densities of the measurement and state, respectively.
The optimal estimation method finds the state that maximizes where S e is the measurement uncertainty covariance matrix of the input brightness temperatures.A priori information are included by x a as the a priori state, that is, the most likely state prior to the measurements, for example, the mean of a climatology, and by S a as the covariance matrix of the a priori values.F is the forward model introduced in Equation 1.The maximum a posteriori solution is given by  x as (Rodgers, 2000) x = Ŝ( with the associated uncertainty where K is the Jacobian at  x of the forward model F. For Equation 6 to be fully valid, we have assumed that the forward model is only moderately non-linear and that the parameter distributions are Gaussian.This is a valid assumption for most of the parameters even though, e.g., LWP and SIC may differ in their distributions from a Gaussian on large scale which could be exploited in future works.

Iterative Procedure
We implemented automatic differentiation (Revels et al., 2016), which allows efficient computing of Jacobians.Then, the cost function Equation 5 is minimized using the Levenberg-Marquard algorithm where the step of the iteration from x i to x i+1 is controlled by the parameter γ i that is either reduced or increased with each step depending on the difference of the cost functions (Rodgers, 2000): We further constrain our parameter space to physically meaningful values (e.g., only positive values of TWV or SIC).If, during the iteration, such nonphysical state vectors are obtained, we re-start using a different startguess (e.g., small but positive values of TWV or SIC).To define a convergence criterion we follow Rodgers (2000) and test whether   2 = ( − +1) ⊺ Ŝ−1 ( − +1) is less than the number of retrieval parameters.

A Priori State x a
As a priori state, we use climatological monthly data based on different data sets, all on a regular lat-lon grid of 0.25°.In the retrieval, the values from the closest grid cell are chosen and linear interpolations between subsequent months are performed.For LWP, WSP, SST, SIC and TWV we use monthly mean values based on ERA5 reanalysis data from 1990 to 2019.From that data, we also obtain 2-m temperature T2m that is then converted to T sa by a seasonally varying regression given in Nielsen-Englyst et al. ( 2021).Snow depth (SND) a priori estimates are based on the merged Warren-AMSR2 climatology (Section 2.4).For multiyear ice fraction we use the multiyear ice satellite product based on microwave radiometer and scatterometer data (Section 2.3) and compute monthly averages for the years 2013-2022.The a priori snow-ice interface temperature of FYI is based on monthly data of vertically polarized TB at 6.9 GHz (TB6V) from 2013 to 2022 of daily AMSR2 data, which is modified for multiyear ice according to an empirical fit (Equation 16in Tonboe et al. (2011)).Using the a priori multiyear ice fraction, a weighted average of T si is calculated.Due to different time spans (30 years of ERA5, but 10 years of TB6V), inconsistencies might arise.Thus, wherever T si is below T sa , we assume an isothermal snowpack with T si = T sa as a priori estimate.

A Priori Covariance S a
For the time being, we assume no correlations yielding a diagonal covariance matrix.The values that we use are given in Table 2 in the last column.The values are based on Scarlat et al. (2017Scarlat et al. ( , 2020) ) and are chosen large enough for the new retrieval parameters to allow deviations from the a priori state.

Start Guess
The startguess x 0 is given by the a priori state vector, except for SIC, Multi-year ice fraction (MYIF) and T si .
The two ice parameters are computed using the NASA TEAM (Cavalieri et al., 1984(Cavalieri et al., , 1997) ) algorithm with tie points originally developed for SSMI sensors (Comiso et al., 1997) using the brightness temperatures at 36.5 and 18.7 GHz.For T si we apply the fit given in Equation 16in Tonboe et al. (2011) using the vertically polarized brightness temperature at 6.9 GHz as input.

Effective Measurement Uncertainty S e
The covariance matrix S e contains measurement but also model uncertainties.Following Maahn et al. (2020), we denote the forward model uncertainty as S b and use an effective measurement uncertainty S e = S y + S b , where S y accounts for the radiometric noise (NEΔT).We assume 1 K for each frequency and polarization for S y and no correlations.In order to get an estimation of the variability lost by fixing the model parameters, we perform a sensitivity analysis to obtain S e .The effect on TB of certain parameters depends also on the retrieval parameters.Thus we scan the retrieval parameter space.For example, changes in ice salinity result in larger changes of computed TB if the ice temperatures are high.Therefore, we vary T sa from 238 to 270 K in steps of 1 K and SND from 0 to 40 cm in steps of 1 cm.The snow-ice interface temperature is chosen from a normal distribution with 1 K standard deviation around the value one obtains assuming thermal equilibrium (using ice and snow conductivities of 2.1 and 0.31 W K −1 m −1 , respectively).The analysis is done separately for first-year and multiyear ice.
We run 1,000 simulations for each tuple of temperature, SND and ice type, in which almost all layer parameters are randomly chosen out of normal distributions that are truncated to constrain the parameter space to realistic values.The standard deviation σ as well as upper and lower bound of these distributions can be found in Table 3.The mean values are the ones given in Table 1.We here focus on the surface emission model.Therefore, for estimating surface model uncertainties, downwelling brightness temperatures reflected at the surface are fixed to typical Arctic values for the six frequencies from 6.9 to 89 GHz (T D = (4, 5, 12, 22.5, 30, 70)K).In the actual   The standard deviation of the resulting brightness temperatures is taken as a measure for the uncertainty that is introduced by fixing the parameter as described above.
With this method we obtain a diagonal covariance matrix S b .Note that, in principle, the uncertainties can be correlated, resulting in nonzero off-diagonal terms in S e , which would affect the weighting of the different TB channels in the cost function.Estimating these correlations is not straightforward, amongst others, a correlation may be conditioned on assumed values for other parameters, as discussed in, for example, Cimini et al. (2018).For the time being, we therefore set the off-diagonal terms to zero.

Additional Sources of Uncertainty
The simulations serve both the understanding and identification of cases with higher uncertainty as well as giving a quantitative measure to adapt the effective error covariance matrix so that it reflects the uncertainties of the model more accurately.However, it is important to keep in mind that not all restrictions of the model were simulated in this approach.There are additional sources of uncertainty.The most notable are: • The number of layers (and thus the possibilities to resolve gradients of, e.g., temperature).

•
The scattering mechanism, that is the parameterization of the scattering coefficient, as well as scattering anisotropy.

•
The mixture model to calculate the effective permittivity of a composite material.
• Snow wetness which is important if temperatures are high (we consider freezing conditions only).
• Roughness of interfaces which we expect to have the largest effect on the higher frequency at 89 GHz.
• Uncertainties can arise for the atmospheric microwave emission model, which was not specifically developed for the polar atmosphere.Thus there are unknown model uncertainties.In addition, one uncertainty is the cloud temperature that is now estimated to be half-way between the surface and freezing point temperature.From a simulation run of the model using the MOSAiC time series as input with cloud temperatures one time fixed at −30°C and one time at 10°C we conclude that the low frequencies have almost no sensitivity to the cloud temperature with mean TB differences of a few millikelvin.We observe the largest effect for 89 GHz which is on average still below 1 K.

Forward Model Evaluation Using MOSAiC In Situ Data
In order to evaluate the forward model, we use MOSAiC observations from October to April as "ground truth" input to the model and compare the simulated brightness temperatures with satellite observations.Additionally, we simulate the uncertainty of the ground truth by running the model for each TB measurement 100 times, where the ground truth values are taken from respective normal distributions with standard deviations (representing spatial and/or temporal variability) given in Table 4. Figure 1 shows the time series of measured and modeled brightness temperatures at the lowest and highest frequency (6.9 and 89 GHz) as well as their probability density distributions and differences.For the other frequencies the reader is referred to Figure B1 in Appendix B. The standard deviations of the 100 simulations are shown as shaded areas around the mean modeled TB in Figure 1.The temporal evolution is well represented in the modeled data.We note less agreement in the horizontal polarization, where ΔTB shows higher values and also the squared Pearson coefficient as a measure for the correlation (r 2 ) is lower (e.g., 0.92 compared to 0.61 for 6.9 GHz).As the satellite measures close to the Brewster angle, the surface reflection (mainly dependent on the difference of the refractive indices of surface and air) is stronger for horizontally polarized TBs and this polarization is therefore more sensitive to changes in snow and ice properties which are not fully represented in our model.Thus, we expected less agreement for horizontal polarization.
In both polarizations, the differences between measurements and model are larger for the higher frequencies.The excellent agreement for vertically polarized TB at 6.9 GHz is due to its high correlation with T si .We also note changes over the season in the differences of modeled and measured TBs, especially at 89 GHz where ΔTB is, for example, lower in March/April compared to November/December.One possible explanation is snow metamorphism: 89 GHz is very sensitive to the upper snow layer microstructure,  suggesting that our model setup is more representative of a snowpack in the late winter season than at the beginning of winter.
Similar correlations between modeled and measured TB are observed when using the RRDP data set (Pedersen et al., 2019) (using Operation Ice Bridge SND and ERA5 data for the other variables as ground truth, see Appendix D).Due to the overall good agreement between simulations and observation, and because the observed biases show a seasonal dependence, we do not make a constant bias correction of the forward model.Some of the deviations from measurements, however, cannot be explained by the input uncertainties and are related to model (parameter) uncertainties.Here, more insight will be gained in the next Section.

Uncertainty Analysis: Covariance Matrices
As described in Section 3.2.5 a sensitivity analysis can be used for evaluating the non-linearity of the model and to estimate the model uncertainties introduced by fixing the model parameters.We use the standard deviation of all simulations (for all tuples of ice type, temperature, and SND) as a measure of the model uncertainty.The standard deviations are plotted as violin plots in Figure 2 separated by ice type, frequency and polarization.For each frequency, the median standard deviation (over both ice types) is listed in Table 5 and is used in the effective measurement covariance matrix S e in the inversion method to include the model uncertainties as described in Figure 1.Modeled and measured brightness temperatures at 6.9 and 89 GHz both vertically polarized (left column) and horizontally polarized (right column).First row shows satellite-based measurements of brightness temperature co-located to Polarstern (dashed lines) and modeled data using ground truth measurements as input (solid lines).The modeled data is given with a 1−σ uncertainty estimate derived from the ground truth uncertainties, see Table 4.In the center row, the distributions of modeled and measured brightness temperatures are plotted as probability density plots.Squared correlation coefficient (r 2 ) and mean absolute error (mae) are shown as annotations in the Figure .The last row displays the difference between measured and modeled data, including uncertainty ribbon.
Section 3.2.5.We can see that for the low frequencies the uncertainty of modeled horizontally polarized TB is higher than for vertically polarized TB while it is the other way around for 36.5 and 89 GHz.For the low frequencies, this was expected due to the Brewster effect.
In general, the uncertainty at vertically polarized TB increases with frequency and is as high as 7 K for 89 GHz.Interestingly, the range of standard deviations is higher for vertically polarized TBs.Depending on the frequency and polarization, different parameters are the dominant sources of uncertainty.For example, ice thickness and corresponding uncertainties of T si , especially for large SND and high T sa are important for vertical polarization and low frequencies, while snow density plays a major role for low-frequency horizontally polarized TB.High frequencies like 89 GHz are especially sensitive to the snow correlation length.
The interplay of parameter uncertainties is captured by our approach, for example, changes in salinity have a larger impact at high temperatures.In order to obtain reasonable results, a sensible choice of simulation ranges is key.
We observe about the same order of magnitude and frequency-dependence of the uncertainties compared to the mean differences of brightness temperatures modeled from ground truth versus satellite measurements (described in the previous Section 4.1).Thus, the deviations between model and measurements can partly be attributed to the parameters of the model setup, that is, better agreement could be obtained by adapting the model parameters accordingly and in a plausible range.We note however that the modeled uncertainties in this section for horizontally polarized TB at 89 GHz are lower than the observed deviations.It is likely that additional model constraints play a role, which were not simulated here (e.g., scattering mechanisms).

MOSAiC Time Series
When we apply the inversion method on the co-located satellite data during MOSAiC, we retrieve parameters as shown in Figure 3.For all cases, the convergence criterion is met.The uncertainty estimate given by Equation 7is shown as the 1-σ shaded area around the maximum a posteriori solution (blue line).It is clear that the method does not reproduce the a priori (shown as green line) in most cases, but instead converges to a different value (most drastically visible for SND).Most of the time, the observed values lie within the retrieval uncertainty.Notable is the sensitivity to LWP, where many Frequency (GHz) Median σ (K) (H/V)  events with high LWP are detected by the retrieval method.The variability of T sa is not captured well in the retrieval, and the values stay close to the a priori.This is likely due to the lower weights on the high frequencies (high values in the effective measurement covariance matrix S e ) that are more sensitive to that parameter.On the other hand, due to the high weight on 6.9 GHz and strong correlation of this channel with T si , the evolution of this quantity is well captured by the retrieval.
For TWV, the agreement is good.Special events like the two warm air intrusions in April 2020 (visible as peaks in the time series), described in Kirbus et al. (2023) and Rückert et al. (2023), are well represented in the retrieved data.Figure 4 shows the good correlation (squared Pearson coefficient: 0.7) and low bias (0.18 kg m −2 ) between satellite retrieved values of TWV and the ship-based measurements.In comparison to the previously used version of the retrieval (using empirical emissivities (Scarlat et al., 2017)), we note a reduction of bias and increase of correlation (see Appendix C).Similar results for TWV are observed in a comparison to radiosonde measurements during the N-ICE campaign 2015 in the sea ice north of Svalbard, see Appendix E.

Arctic-Wide Multi-Parameter Retrieval
We run the retrieval on satellite swath data of a whole day and the retrieval parameters are then resampled (Gaussian resampling using NSIDC's Polar Stereographic Projection (A Guide to NSIDC's Polar Stereographic Projection, n.d.) with a nominal gridded resolution of 12.5 km by 12.5 km) to obtain maps on a pan-Arctic scale.One example for 10 January 2020 is shown in Figure 5, one for 19 April 2020 in Figure 6, the time when warm and moist air from the North Atlantic intruded into the central Arctic (Kirbus et al., 2023).More examples can be found in Appendix F. The convergence rate is above 98%.
The sea ice extent is well captured (compared to ERA5, not shown).The spatial distributions of multiyear ice fraction and SND are plausible, with higher values in the Canadian Arctic.The distribution of snow-air interface temperatures and snow-ice interface temperatures are within expected ranges.Higher values of TWV and LWP over the open sea are observed as we can expect it.In some areas, the artificial gradient of TWV and LWP with respect to the ice type are still visible, but are a lot less pronounced than for the previously used version of the retrieval (see Figure C1).The warm and moist air intrusion on 19 April 2020 is also well seen in the retrieved elongated filament of high values of LWP and TWV.Some areas over FYI (e.g., Chukchi Sea around Bering Strait, Laptev Sea) exhibit higher values than the reanalysis ERA5 (not shown).
When we investigate the retrieval uncertainty with respect to water vapor (Figure 7) we see clearly higher uncertainties over the sea ice (middle).Compared, however, to the a priori uncertainty of 4.69 kg m −2 we see an improvement also over ice, seen as values smaller than one in the ratio of retrieved to a priori uncertainty (Figure 7, right).Such tests help in the identification of areas of higher uncertainties, in this example the Russian Arctic with younger ice around 150°E.
An additional analysis using the RRDP with confirmed scenes of 100% SIC from October to April (including observations from the predecessor AMSR-E sensor) shows that our retrieval method has an accuracy and precision similar to or even better than many other common SIC algorithms (Ivanova et al., 2015): we retrieve 100% SIC with a mean bias of 0.5% and a standard deviation of 1.5% (note, however, that our method truncates SIC to 100% and does not allow SIC values >100%).
In summary, these maps and the additional ones shown in the Appendix serve as a first check on the retrieval.More research is underway to evaluate the retrieval further, using more reference data sets and investigating the seasonal evolution of the different parameters.The evaluation with MOSAiC observations in Section 4.3.1,however, already demonstrates the good performance of our new multi-parameter retrieval especially compared to the previous version.

Conclusions and Outlook
A combined atmosphere-snow-ice-ocean emission model together with optimal estimation has been used in a multi-parameter retrieval of snow, ice and atmospheric parameters in the Arctic from satellite AMSR2 microwave radiometer observations.Instead of correcting the satellite measurements for weather influences, we make use of the atmospheric information by using an integrated retrieval.Using literature values and in situ data, a forward model has been designed to model top-of-atmosphere brightness temperatures from nine geophysical parameters.
The surface part of the model requires four input variables, namely snow-air interface temperature T sa , snow-ice interface temperature T si , the ice type (either FYI or multiyear ice [MYI]) and SND.Data from the MOSAiC expedition provide valuable descriptions of the "ground truth" that can be used as model input.Despite simplifications regarding model parameters, we can simulate realistic values of brightness temperatures, and their variabilities, resembling satellite observations by using MOSAiC input data.Similar results are obtained using the RRDP with reference sea ice observations.To some extent, differences between forward model and observation can be explained by model parameter choices as we show in the simulation of the uncertainties of the model (Section 4.2).In addition, in the MOSAiC case, we observe a seasonal dependence of the bias with better agreement toward spring (see e.g., Figure 1).More research is under way to attribute the seasonal performance dependence to certain model parameters which can then be refined in the future.
The model is invertible using an optimal estimation method, which provides us with a physically consistent set of nine geophysical parameters and their uncertainties.Although the degrees of freedom in the measurements are less than the number of retrieval parameters, they are mostly in agreement compared to MOSAiC observations.The model inversion by optimal estimation comes at acceptable computational costs, which enables daily retrievals using all available satellite observations from AMSR2 (and its predecessor AMSR-E).Spatial patterns (Figures 5 and F1) as well as the seasonal development (Figure 3) of the retrieved parameters are plausible.This is now also true for the atmospheric parameters TWV and LWP in comparison to the previous model that used constant, empirical emissivities for sea ice with snow.
More comprehensive evaluations of the retrieved parameters can now follow.For example, some areas are persistently featuring high values of LWP over weeks, which is likely due to the surface representation rather than actual atmospheric events.It indicates that the parameterization of FYI in our model setup might need further improvement.Here, more observational data, especially of snow microstructure, could help to improve the model setup.Also, reference data sets from reanalysis and other satellite retrievals are needed.Available and suitable ground observations unfortunately remain sparse.Another challenge is the retrieval under summer melting conditions for which our sea ice forward model would have to be adapted.At our microwave frequencies water in melt

Appendix D: Round Robin Data Package D1. Forward Model Evaluation
We follow a similar approach as with the MOSAiC data sets for the Round Robin Data Package (RRDP, Pedersen et al., 2019), and assume the RRDP data to be a ground truth.Most parameters included in RRDP are used as given therein (SIC from the percentage of open water, snow depth [SND], TWV, LWP, WSP).We set SST to 273.15 K and make an assumption about MYIF, which is set to one if the sea ice thickness variable in RRDP is larger 1.5 m, else it is set to 0. We calculate T sa from T2m in RRDP, and T si is calculated from TB6V with the corresponding fit over MYIF (see Section 3.2.2).
The results of the comparison between modeled TBs and satellite TBs are summarized in Table D1.Similar to the comparison using MOSAiC in situ data, we find high correlations between modeled and satellite TBs, which is decreasing with frequency.Mean absolute error and bias are also again smaller for vertical polarization than for horizontal.Note.The simulated TB's are compared to the corresponding satellite measured TB's included in the RRDP and mean absolute error, bias and r 2 are listed in this Table .The total number of data points is 530 (for each frequency).

Table D1
From the Round Robin Data Package, Ground Truth Data Is Taken as Input for the Forward Model

Appendix F: Arctic-Wide Retrieval: Daily Maps
In order to gain a first look at the seasonal evolution of the parameters and their spatial distributions, we show daily maps of all parameters from 1 day in October, December, and February, see Figures F1-F3.

Data Availability Statement
The AMSR2 Level 1B and Level 1R satellite data in the study is available via Japan Aerospace Exploration Agency's G-Portal https://gportal.jaxa.jp/gpr/?lang=en.The ERA 5 reanalysis data used for calculating a priori climatologies in the study are available at the Copernicus Climate Change Service (C3S) (2023) (Hersbach et al., 2023).The multi-year ice fraction data used for calculating a priori climatologies in the study are available from the University of Bremen's research group https://www.seaice.uni-bremen.de/start/via https://seaice.uni-bremen.de/data/MultiYearIce/.The merged Warren-AMSR snow depth climatology was obtained from the Alfred Wegener Institute Helmholtz Center for Polar and Marine Research (AWI) via ftp://ftp.awi.de/sea_ice/auxiliary/snow_on_sea_ice/w99_amsr2_merged/.The campaign data sets as part of the international Multidisciplinary drifting Observatory for the Study of the Arctic Climate (MOSAiC) are obtained from the following sources: • Total water vapor and liquid water path are available at PANGAEA via https://doi.org/10.1594/PANGAEA.941389(Ebell et al., 2022).• Snow buoy snow depth data is available at PANGAEA via https://doi.pangaea.de/10.1594/PANGAEA.938244(Lei et al., 2021a).Nixdorf et al. (2021).We gratefully acknowledge the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Transregional Collaborative Research Centre TRR-172 "ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC) 3 " (Grant 268020496), the MOSAiC-microwaveRS project (Grant 420499875) and ESA's CIMR-DEVALGO project (4000137493/22/NL/AD).R. T. Tonboe was supported by the ESA project: Performance Evaluation of Arctic Weather Satellite Data (4000136511/21/NL/ IA).Special thanks to Wenkai Guo for providing the classified TerraSAR-X sea ice scenes and JAXA for providing the AMSR2 data.We wish to thank the two reviewers and the editor for their helpful comments.Open Access funding enabled and organized by Projekt DEAL.

Figure 2 .
Figure 2. Sensitivity analysis for effective measurement covariance estimate.Shown are simulated standard deviations of 1,000 simulations for each temperature and snow depth combination, separated by frequency, polarization, and ice type.The violin plots show the 5th-95th percentile.

Figure 3 .
Figure 3. Multidisciplinary drifting Observatory for the Study of Arctic Climate timeseries of retrieved parameters (blue) with their uncertainties (shaded areas), additionally the ground truth observations as described in Table4(gray line) and the a priori data (green line) is shown.Note that we do not display the full observational data sets, but only the part co-located with satellite measurements.The following seven quantities are displayed: SIC: sea ice concentration; MYIF: multiyear ice fraction; T sa : snow-air interface temperature; T si : snow-ice interface temperature; SND: snow depth; TWV: total water vapor; and LWP: liquid water path.

Figure 4 .
Figure 4. Scatter plot of retrieved total water vapor using the emission model versus Multidisciplinary drifting Observatory for the Study of Arctic Climate measurements.Included in the plot are bias, mean absolute error (mae), and the coefficient of determination given by the square of the Pearson correlation coefficient, that is, r 2 .

Figure 5 .
Figure 5. Retrieved parameters, gridded and resampled from daily AMSR2 swath data from 10 January 2020.Shown are sea ice concentration (SIC), multiyear ice concentration obtained by multiplying the multiyear ice fraction with SIC, total water vapor and liquid water path in the first column.The second column displays Arctic-wide maps of the surface temperature as average (weighted by SIC) from the snow-air interface temperature T sa over ice and the sea surface temperature over open water, the snow-ice interface temperature (T si , shown only in areas where SIC > 15%), as well as wind speed (shown only in areas where SIC < 15%) and snow depth (shown only in areas where SIC > 15%).

Figure 6 .
Figure 6.Retrieved parameters, gridded and resampled from daily swath data as in Figure 5 but from 19 April 2020.

Figure C1 .
Figure C1.Map of retrieved liquid water path (daily gridded data) using the model with empirical emissivities for 10 January 2020.The 0.5 contour line (yellow) of retrieved multi-year ice fraction is shown in addition.

Figure C2 .
Figure C2.Scatter plot of retrieved total water vapor using the model with empirical emissivities versus Multidisciplinary drifting Observatory for the Study of Arctic Climate measurements.Included in the plot are bias, mean absolute error (mae), and the coefficient of determination given by the square of the Pearson correlation coefficient r 2 .

Figure D1 .
Figure D1.Density plot of snow depth data.Shown are observations from operation ice bridge included in the Round Robin Data Package (gray), corresponding retrieved values (blue) and the a priori values used in the retrieval (green).

Figure F1 .
Figure F1.Retrieved parameters, gridded and resampled from daily swath data as in Figure 5 but from 10 October 2019.

Figure F2 .
Figure F2.Retrieved parameters, gridded and resampled from daily swath data as in Figure 5 but from 10 December 2019.

Figure F3 .
Figure F3.Retrieved parameters, gridded and resampled from daily swath data as in Figure 5 but from 10 February 2020.

Table 2
Source of A Priori State x a (From Climatology) and A Priori Covariances Given by Standard Deviation σ for All Retrieval Parameters

Table 3 Table
With Parameter Ranges Used in the Simulations to Estimate Model Uncertainty 10 of 29 forward model and retrieval procedure downwelling brightness temperatures are calculated by the atmospheric model.

Table 4 "
Ground Truth" Data Sets Used in the Model Evaluation and Assumed Variability for Each Parameter