© V. Slemzin et al., Published by EDP Sciences 2016

1. Introduction

The influence of the Earth’s atmosphere is very important for spacecraft and scientific instruments operating in near-Earth space. The lifetime of a satellite at altitudes of several hundred km strongly depends on the density of the atmosphere and its variation with solar activity, so the prediction accuracy of atmospheric conditions has a significant economic aspect. The solar short-wavelength (X-ray and EUV) radiation is appreciably attenuated in the upper atmosphere at altitudes below 500 km and fully absorbed below 100 km, which should be taken into account in planning solar observations from satellites in the low-Earth orbits (LEO). On the other hand, due to the high sensitivity of solar observations to atmospheric absorption, the data obtained in the presence of absorption can be used for validation and improvement of atmosphere models.

Spectrometric investigations of the transparency of the Earth atmosphere in short wavelengths as a function of altitude began with rocket launches and first satellite flights during the 1950s to 1970s with the aim of studying absorption of the atmosphere under different observational conditions (see the review by Smith & Hunten 1990 and references therein). For example, Hinteregger (1962) describes such measurements at low zenith angles in the spectral band 28.3–121.6 nm for altitudes of 100–200 km. Later, measurements were carried out from satellites during occultations at the sunset or sunrise parts of the orbit. Below the 100 km level, the measurements of radiation flux variations during solar or stellar occultations are used for monitoring of ozone, nitric, and halogen molecules as well as other minor atmosphere constituents in many spectral bands from UV to infrared (experiments Nimbus-7/TOMS, Bowman & Krueger 1985; SAGE I–III, Chu et al. 2003; Damadeo et al. 2013; UARS/HALOE, Russell et al. 1993; ENVISAT/GOMOS, Bertaux et al. 2010). The solar occultation technique was widely used for remote sensing of the atmospheres of planets and satellites such as Mars, Saturn, and Titan (Smith et al. 1983; Smith & Hunten 1990; Krasnopolsky et al. 1991; Korablev 2002; Maltagliati et al. 2013). The foundation of the occultation method including applications to radio waves and the GPS-based localization were treated in detail in the reviews of Elliott (1979), Smith & Hunten (1990) and in the book “Occultations for Probing Atmosphere and Climate” (Kirchengast et al. 2004).

Systematic measurements of the absorption of solar EUV radiation in the Earth’s atmosphere are not numerous. Since 2002, the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite is monitoring the EUV solar radiation. The Solar EUV Experiment (SEE) on board TIMED (Woods et al. 2005) provides occultation measurements of solar irradiance through the upper atmosphere at different EUV wavelengths, but due to telemetry constraints measurements have been performed at only one single altitude on each orbit. Thus, the whole profile in the range 100–400 km is summarized from measurements on many orbits.

Due to the integration of the measured absorption along the line-of-sight (LOS), occultation measurements have the evident drawback, as a method for validation of atmosphere models, of not permitting to obtain local densities directly. At the same time, this method provides self-calibration of the data, high height resolution, and short time of scanning through a wide range of altitudes which allows measuring the atmosphere absorption on a wide spatial scale.

The ordinary method of solar occultation uses measurements of the total flux integrated over the full solar disk as a function of the so-called tangential altitude h (see Fig. 1), i.e. the nearest distance between the LOS and the Earth surface (sub-LOS point). Therefore, the transmittance profile is averaged over the altitude interval corresponding to the angular size of the Sun (~0.53°), which for a typical satellite orbit with an altitude of about 600 km varies from 1–2 km at h ~ 400 km to ~20 km at h ~ 200 km. The use of high-resolution solar images, instead of integrated fluxes, permits to diminish this uncertainty to the size of one pixel, which corresponds to only several arc seconds, thus increasing the spatial resolution to scales below 1 km. Additionally, bright active regions, where the flux is the most variable, may easily be isolated by intensity and excluded from the integration (Slemzin et al. 2003).

Figure 1. Scheme of the occultation measurements. R0 – radius of the Earth, Horb – altitude of the satellite orbit. Other notations are given in the text.

Since the launch of the first spaceborne EUV solar telescope Phobos/TEREK in 1988 (Zhitnik et al. 1998) and the next CORONAS-I/TEREK-C telescope in 1994 (Sobelman et al. 1996), it was realized that atmospheric absorption effects are very important for regular solar imaging in the EUV spectral range. One can completely avoid the absorption effects if the spacecraft is in the Lagrangian L1 point like the SOHO observatory, or sufficiently diminish them in time if the satellite operates in a sun-synchronous orbit. Due to precession, all low-Earth orbits (LEO) have seasonal periods, when the line-of-sight crosses the Earth atmosphere. In most cases, observational programs of solar instruments do not include such periods in order to avoid distortion of the results by attenuated data. However, the scientific programs of the CORONAS-F/SPIRIT (Slemzin et al. 2005), CORONAS-Photon/TESIS (Kuzin et al. 2009), and PROBA2/SWAP (Seaton et al. 2013) telescopes foresaw special periods of observations of the Sun in occultations to investigate the response of the upper atmosphere to long-term or short-term variations of solar and geomagnetic activity. The first results of determination of atmospheric transmittance carried out by SPIRIT in the 0.842, 17.5, and 30.4 nm spectral bands as a function of tangential altitude in solar maximum (October–November 2001) were published by Zhitnik et al. (2003).

So far, several models intended for the calculation of characteristics of the upper atmosphere such as densities of the main constituents and temperature have been developed. These models allow calculating densities of the main atmosphere constituents at different altitudes, for given geographic coordinates, date, time, and solar and geomagnetic conditions. The first models (e.g. Jacchia 1971) were based on satellite drag data only. The series of the Mass-Spectrometer-Incoherent-Scatter (MSIS) models progressively improved accuracy through the use of the mass spectrometer and incoherent scatter radar data. These models are named MSISE-83, MSISE-86, MSISE-90 (Hedin 1983, 1987, 1991) and the latest NRLMSISE-00 (Picone et al. 2002), which allows calculating densities of the main constituents and kinetic temperatures at altitudes from 0 to 1000 km. The NRLMSISE-00 model can be run online at http://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php.

Comparisons of the NRLMSISE-00 model calculations with the drag measurements on the satellites ARGOS in the altitude range 150–200 km (2000–2001, Aksnes et al. 2007) and CHAMP at about 400 km (2002, Zhou et al. 2009) have shown that under quiet geomagnetic conditions the temperature and density height distributions given by the model were in satisfactory agreement with the measurements. However, it was found, that the model was unable to describe well the rapid enhancements in the thermospheric mass density and temperature during geomagnetic storms in October–November 2003. Density enhancements during storms are very important for the LEO spacecraft such as the International Space Station, because they can seriously alter the orbit and reduce altitude and useful lifetime of the mission.

A series of semiempirical Drag Temperature Models (DTMs) were developed since 1978 (Barlier et al. 1978). The latest advanced model DTM-2013 has been developed under the ATMOP FP7 project with the goal of meeting the requirements for operational orbit computations (Dudok de Wit & Bruinsma 2011; Bruinsma et al. 2012; Bruinsma 2015). This model was constructed using drag observations of 10 satellites, gathered during the period of 1969–2011. The data include drag measurements obtained with accelerometers on board the CHAMP (Bruinsma et al. 2004), GRACE (Bruinsma & Forbes 2008), and GOCE (Bruinsma et al. 2014) satellites and use the F30 (30 cm radio flux) solar activity proxy which proved to be more representative than the F10.7 index (Dudok de Wit et al. 2014). A detailed description of the DTM-2013 model and its comparison with the previous models DTM2009 (Bruinsma et al. 2012) and JB2008 are given in Bruinsma (2015). The online version of the model can be run at http://www.atmop.eu/index.php/models). In our work, we employed the model DTM-2013 for comparison with absorption measurements.

Measurements of atmospheric attenuation of solar EUV radiation constitute an independent way of validating these atmosphere models at different levels of solar and geomagnetic activity. The local extinction coefficients can be retrieved mathematically from the measured absorption and compared with the model predictions. Here we describe the results of occultation measurements of the atmospheric transmittance in the 17 nm range at different altitudes in the quiet state of the atmosphere and during the development of a geomagnetic storm. The data have been obtained with the SPIRIT and TESIS telescopes on board the CORONAS satellites as well as with the SWAP telescope on board PROBA2. We compare the transmittance profiles and the retrieved extinction coefficients as a function of altitude with calculations by the NRLMSISE-00 and DTM-2013 atmosphere models.

2. Data and analysis

Figure 1 shows the scheme of solar observations from satellites on LEO orbits in occultation mode. The line-of-sight (LOS) crosses the atmosphere at altitude h (tangential height) and forms an angle z with the direction to the Earth’s center. The rays from the satellite to different points on the solar disk deviate from the central ray by angles θ within half of the angular size of the Sun (~0.53°). This corresponds to variations dh(θ) of the tangential height h.

Due to the finite size of the Sun, the highest and the lowest rays differ in height by Δh ≈ α(R₀ + H_orb)cos z, where α is the angular size of the Sun, H_orb is the altitude of the orbit, z is the angle between LOS and direction to the Earth’s center, and R₀ is the Earth’s radius. For instance, if H_orb = 600 km and h = 300 km, ∆h ≈ 18 km.

Transmittance T of the atmosphere at a specific wavelength λ along the ray corresponding to a LOS height h and angle θ is given by the integral over the LOS$$T\left(\lambda ,h,\theta \right)=\exp \left(-\underset{\mathrm{LOS}\left(\theta \right)}{\overset{}{\int }}{k}_{\lambda }\left(x,\theta \right)\mathrm{d}x\right),$$(1)where k_λ(x) is the wavelength-dependent extinction coefficient. Generally, extinction (or attenuation) of the light is caused by three processes – absorption, refraction, and scattering of light in the medium. Due to the very small density of the atmosphere at altitudes above 100 km (more than seven orders less than that on the ground), refraction and scattering produce negligible effects, so photoabsorption caused by photoionization of the atmospheric constituents through solar EUV radiation dominates. Thus, the extinction coefficient can be expressed in the form:$$k_{\lambda }\left(x,\theta \right)=\sum _i{d}_i\left(h\left(x\right)+\mathrm{dh}\left(\theta ,x\right)\right){\sigma }_i\left(\lambda \right),$$(2)where d_i(h) and σ_i(λ) are the density of the ith constituent (provided by the atmosphere model) and its photoabsorption cross-section at wavelength λ.

In general, using solar imaging telescopes one can obtain transmittance profiles in two ways: by integrating the solar flux over the whole disk (the integrated mode), or by dividing the disk into many horizontal strips corresponding to discrete heights and by constructing the total transmittance height profile from partial profiles for these strips (the imaging mode).

In the integrated mode, which was used in this work, the signal-to-noise ratio is the highest, but the resulting profile is averaged over the height divergence interval. However, it can be shown, that for a quasi-uniformly radiant solar disk (which is the case in the EUV 17 nm band) the resulting profile coincides with that for the central ray within ∆h ~ 0.5 km. In the imaging mode, the height resolution is defined by the size of the resolved image cell which can be set as small as possible only limited by the signal level and by the growing influence of local temporal variations of the solar emission. The tests have shown that by using the imaging mode during quiet solar condition one can easily achieve a geometrical height resolution better than 1 km (Slemzin et al. 2003).

It is worth noting that even in the integrated mode the height divergence due to the solar size is less than the atmosphere scale height H_s. Calculations with DTM-2013 model for a 600 km orbit in the tangential height range h = 200–500 km give values of H_s that vary from 42 to 104 km at solar maximum (the SPIRIT 1st data set) and from 24 to 120 km at minimum (the TESIS data set). Taking into account smooth dependence of density on height in the given range, the finite solar size cannot produce significant errors in the transmittance profiles.

We constructed transmittance profiles using solar images obtained in the sunset or sunrise sections of the orbit. To calculate the positions of the satellites and the sub-LOS points (at the Earth surface) at the moment of image acquisition we used two-line element (TLE) ballistic data provided by the Joint Space Operations Center of the US Strategic Command (https://www.space-track.org), processed with the TRAKSTAR code developed by T.S. Kelso (http://celestrak.com/software/tskelso-sw.asp). As a rule, the TLE catalog contains one or two ballistic data sets per day for a given satellite from which we selected the nearest data set in time to the measured image. For quiet Sun periods (the TESIS and SWAP data sets), it resulted in an accuracy in the determination of the tangential height of ~25–50 m (this value was derived from the difference by two subsequent TLEs for the same instant of time). Even at solar maximum (the SPIRIT case), the error in h did not exceed 100 m.

The initial images were processed by standard procedures for each instrument which included corrections for the background, the flat-field distortion, and the nonlinearity. Then the intensity was integrated over the disk area excluding the corona above the limb where the EUV flux can significantly vary in time due to solar activity. To obtain transmittance profiles, the integrated light curves were determined as a function of tangential height and then normalized by the mean value of intensity above 500 km. At this altitude the absorption is comparable to or less than the relative variation of the solar flux due to solar activity and can therefore be neglected. In the studied cases, the dispersion of the solar flux at h > 500 km in the considered spectral range did not exceed 0.3% at maximum activity (SPIRIT) and 0.1% at minimum activity (TESIS and SWAP).

We validated the models by confronting them with the atmosphere absorption observed at different levels of solar activity (near maximum, at minimum, and in the growing phase of the solar cycle) under quiet and disturbed geomagnetic conditions. In the quiet or slightly disturbed geomagnetic state, we analyzed the following data sets: the SPIRIT images in the 17.5 nm band obtained on March 5, 2002, between 02:45:35 UT and 03:11:48 UT (data set 1 with cadence 32 s near maximum of solar cycle 23, geomagnetic index ap = 32 nT); the TESIS images in the 17.1 nm band (near solar minimum between cycles 23 and 24, July 17, 2009, 13:45:30 UT–July 19, 2009, 15:20:35 UT, cadence 5 min, ap < 3 nT); the SWAP images in the 17.4 nm band (the growing phase of cycle 24, December 22, 2010, 08:27:19 UT–08:34:39 UT, cadence 20 s, ap = 0). The response of the models to disturbed geomagnetic conditions was studied during the development of a geomagnetic storm using the SPIRIT data set 2 obtained before the storm on 23 May 2002 01:54:20–02:20:33 UT (ap = 12 nT, Dst = −1 nT), and data set 3 during the storm, 15:59:59–16:48:56 UT (ap = 236 nT, Dst = −75 nT). The SPIRIT and SWAP data were collected in the sunrise parts of single orbits, the TESIS data were combined from the sunrise and sunset parts of 30 orbits. Spectral functions for all instruments have half-widths of 0.6–1 nm which is wider than the differences in the peak wavelengths, so they operated in similar wavelength bands.

Figure 2A presents the temporal positions of these data superimposed on the evolution curve of solar activity as testified by the sunspot cycle taken from the WDC-SILSO at the Royal Observatory of Belgium in Brussels (http://sidc.oma.be/silso). Figure 2B shows the transmittance profiles for SPIRIT, TESIS, and SWAP corresponding to solar maximum of the cycle 23, and the minimum and growing phase of the solar cycle 24. Points in Figure 2B correspond to single measurements, solid curves to the fitted profiles. With growing solar activity, the increase of the transmittance from ~0 to ~1 occurs at higher altitudes and extends over a wider altitude range. Difference in transmittance is maximal at the heights 220–350 km, where the atmosphere absorbs the major part of solar EUV flux. The map in Figure 3 shows the geographic locations of the sub-satellite and sub-LOS points for the input data. During the measurements of transmittance profiles, the sub-LOS points shift in latitude and longitude by tens of degrees, so the obtained profiles characterize the distribution of atmospheric densities over a wide spatial range.

Figure 2. (A) Sunspot numbers in solar cycles 23 and 24 and observation times of the SPIRIT, TESIS, and SWAP data (explanations in the text); (B) transmittance profiles for SPIRIT (solar maximum, data set 1, red line), TESIS (solar minimum, black line), and SWAP (rising phase of cycle 24, blue line).

Figure 3. Map of occultation locations: red symbols – locations of sub-satellite points for SPIRIT (data set 1), TESIS and SWAP, blue symbols – locations of sub-LOS points.

3. Comparison of experimental transmittance profiles with model calculations

Figure 4 displays the height profiles of extinction coefficients at 17 nm wavelength of the main atmospheric constituents in the height range 200–500 km – atomic oxygen, molecular nitrogen, molecular oxygen, and helium – calculated with the atmosphere models NRLMSISE-00 and DTM-2013 at solar maximum, representative for the SPIRIT data set 1, and at solar minimum, representative for the TESIS data set. The photoabsorption cross-sections for the dominant atomic components (O, N, He) were taken from Verner et al. (1996) and Fennelly & Torr (1992), for O₂ and N₂ from Verner et al. (1996) and Cole & Dexter (1978).

Figure 4. Height profiles of extinction coefficients of the main atmospheric constituents calculated with the NRLMSISE-00 and DTM-2013 models for the EUV wavelength 17 nm in the quiet geomagnetic state. (A) At solar maximum (the SPIRIT data set 1), (B) at solar minimum (the TESIS data set).

A comparison has shown that in the spectral interval 5–35 nm the difference between the data from different sources does not exceed 2% for O and N and 10% for O₂ and N₂ which is close to the theoretical data accuracy (5–8%). Extinction coefficients for different components calculated with densities provided by both models depend differently on height and level of solar activity. The total extinction calculated with the NRLMSISE-00 model at solar maximum is higher than that calculated with the DTM-2013 model (the ratio averaged over the most important height range 200–500 km is 1.49), at solar minimum the former model gives the smaller value (the ratio is 0.74).

The transmittance profiles corresponding to the above-mentioned SPIRIT, TESIS, and SWAP data were calculated according to Eqs. (1) and (2). The transmittance along LOS as a function of height was determined by numerical integration of the extinction under the assumption of spherical symmetry of the atmosphere. In Figure 5, the modeled profiles for the SPIRIT, TESIS, and SWAP data sets are compared with those derived from measurements.

Figure 5. Comparison of the transmittance height profiles, presented in Figure 2B, with those calculated by the NRLMSISE-00 and DTM-2013 models. The black lines display the averaged measurement points (for TESIS individual points are also shown), the red/blue lines show simulations by the NRLMSISE-00 and DTM-2013 models. In the SPIRIT and SWAP cases, the error bars correspond to the transmittance determination errors, in the TESIS case to the dispersion of data averaged over multiple profiles.

The difference between the measured and modeled profiles can be expressed in mean displacement of the transmittance height profile (in km) averaged over the transmittance range 0.1–0.9 (Table 1). The profiles calculated with the NRLMSISE-00 model are closer to the data fitting functions for periods of high solar activity (SPIRIT set 1) and moderate activity (SWAP), the DTM-2013 model profile is closer for solar minimum (TESIS), although the mean displacements in the SPIRIT and TESIS cases are below the data errors, which included the intensity measuring error and variation of the unabsorbed solar flux during measurements.

4. Comparison of the retrieved extinction coefficients with the modeled values

As it was shown in Section 2, a comparison of the modeled transmittance height profiles with those derived from measurements gives an indication of how the models describe the large-scale distribution of atmospheric densities, but does not permit a detailed estimation of the accuracy of their local modeled values. For this purpose, in the cases of SPIRIT and SWAP we calculated the local extinction coefficients as a function of height at sub-LOS points with both models and compared them with the values retrieved from the measured profiles.

We have developed an approach for solving the inverse problem of retrieval of the extinction coefficients from measured profiles based on Tikhonov’s theory of ill-posed problems and regularization method (Tikhonov & Arsenin 1977), which allows overcoming the inevitable limitations inherent in the foregoing model-based consideration; in particular, this approach makes it possible to retrieve variations of the absorption coefficient with fine altitude resolution. The basic principles of this method and its application to the analysis of the attenuation profiles are presented in Annex A. While the attenuation profiles characterize the density distribution in the atmosphere in the wide scale, the extinction coefficients are directly defined by the local density near the sub-LOS point at a given height, so they are very useful for validation of the atmosphere models. Figures 6 and 7 illustrate a comparison of the extinction coefficients (at the sub-LOS points) retrieved by the method described in Annex A with those calculated using the densities provided by the NRLMSISE-00 and DTM-2013 models.

Figure 6. Comparison of the modeled and retrieved extinction coefficients as a function of height for the SPIRIT data set 1 (March 5, 2002, left) and the SWAP data set (December 22, 2010, right). Black line – the data retrieved from measurements, red/blue lines – calculations by the NRLMSISE-00 and DTM-2013 models.

Figure 7. Variation of the retrieved and modeled extinction coefficients during the development of the geomagnetic storm of 23 May 2002: (A) before the storm (the SPIRIT data set 2, Dst = −1 nT), (B) 15 h later (the SPIRIT data set 3, Dst = −75 nT). The line colors are the same as in Figure 6.

4.1. Quiet atmosphere at different levels of solar activity

Figure 6 shows a comparison of the extinction coefficients as a function of height retrieved from the SPIRIT data (set 1, March 5, 2002) and the SWAP data (December 22, 2010) with those modeled with NRLMSISE-00 and DTM-2013. The SPIRIT data correspond to a slightly disturbed geomagnetic state (the geomagnetic index ap = 32 nT), the SWAP data to the quiet state (ap = 0 nT).

As seen in Figure 6, at altitudes from 220 to 460 km, in the SPIRIT case the NRLMSISE-00 model underestimates the extinction coefficients of ~30%, whereas the DTM-2013 model agrees with the data within 5%. In the SWAP case, the models show mean deviations from the data that reach 48% (NRLMSISE-00) and 39% (DTM-2013) at heights above 300 km. It can probably be explained by the fact that the PROBA 2 satellite operates in a dawn-dusk orbit, so the sub-LOS points for SWAP are located near the terminator where density gradients and the uncertainty of the simulation are the largest. The accuracy of the retrieved extinction coefficients in all cases was about 5% due to the self-consistent mathematical procedure (see Annex A2), which is better than the accuracy of the individual data values.

4.2. Disturbed atmosphere (geomagnetic storm)

Figure 7 shows the variation of the retrieved extinction coefficients during the development of the moderate geomagnetic storm on May 23, 2002, 01–15 h UT (the SPIRIT data sets 2 and 3). During this period, the ap index rose from 12 nT to 236 nT, and the Dst index fell from −1 nT to −75 nT. The storm began at ~9 h UT (http://omniweb.gsfc.nasa.gov). The largest variation of the extinction occurred at the heights 300–450 km (increase at 400 km by a factor of 1.7) which was probably produced by density enhancements in the ionosphere associated with magnetic disturbances during the storm (Hedin et al. 1977).

The NRLMSISE-00 model describes well (within the data accuracy) the variation of the extinction coefficient in the quiet state before the storm, but fails to describe its growth during the storm (difference about 70%). The DTM-2013 model describes well the behavior of the extinction coefficients above 300 km height in the quiet state (accuracy ~12%) and under storm conditions (~1%). Below 300 km the extinction coefficient given by this model exceeds the one retrieved from measurements, at 200 km by a factor of 2. At the moment, the reasons for such a discrepancy are not clear.

5. Summary and conclusions

1. The study of the atmospheric attenuation of solar EUV radiation using telescopes and photometers on board LEO satellites during occultations provides valuable information about the properties and condition of the Earth’s upper atmosphere and can be used for validating atmosphere models. The EUV wavelength band is well suited for the investigation of the atmosphere at heights 200–500 km where the influence of solar and geomagnetic activity is the largest.

2. The transmittance height profiles during solar occultations obtained from the total integrated flux over the whole disk have a limited height resolution of 10–20 km due to averaging over the finite angular size of the Sun. However, under the assumption of a quasi-uniform distribution of brightness over the solar disk (which is valid in the EUV 17 nm band), the resulting profile coincides with that of the ray from the solar center within ±0.5 km. Earlier (Slemzin et al. 2003) it was shown that the height resolution of the profiles can be increased to less than 1 km if using solar images segmented into narrow strips.

3. For the analysis of the atmospheric absorption we used solar observation data obtained during orbital occultations by the telescopes CORONAS-F/SPIRIT (17.5 nm), CORONAS-Photon/TESIS (17.1 nm), and PROBA2/SWAP (17.4 nm). The selected data were obtained at different levels of solar activity: at the maximum of cycle 23 (SPIRIT, March–May 2002), in the deep minimum between cycles 23 and 24 (TESIS, July 2009), and at moderate solar activity during the growing phase of the cycle 24 (SWAP, December 2010).

4. We compared the obtained transmittance profiles at the heights 200–500 km with simulations based on the atmosphere models NRLMSISE-00 and DTM2013. The profiles calculated with the NRLMSISE-00 model are closer to the data fitting functions for the periods of high solar activity (SPIRIT set 1, the mean displacement about 1 km) and during a moderately active Sun (SWAP data set, displacement 4.5 km), the DTM-2013 model for solar minimum (TESIS data set, displacement 1.7 km).

5. Extinction coefficients were retrieved with an accuracy of about 5% from the measured transmittance profiles by solving the inverse problem using Tikhonov’s theory of ill-posed problems and regularization method based on the generalized discrepancy principle. It was found that at altitudes between 220 and 460 km, under slightly disturbed magnetospheric conditions (in the SPIRIT case) the NRLMSISE-00 model underestimates the extinction coefficients by ~30%, whereas the DTM-2013 model agrees with them within 5%. In the SWAP case (quiet magnetosphere), the mean deviations from the data at heights above 300 km are higher (48% for NRLMSISE-00 and 39% for DTM-2013), which can probably be explained by the uncertainty of the simulations in regions of high density gradients.

6. The response of the upper atmosphere to geomagnetic disturbances was demonstrated by the comparison of the retrieved extinction coefficients for the SPIRIT orbits 4508 and 4524 (23 May 2002, with starting times 01:54 and 15:59 UT) corresponding to transition from the quiet magnetospheric conditions (ap = 12 nT, Dst = −1 nT) to a moderate storm (ap = 236 nT, Dst = −75 nT). Under geomagnetic storm conditions we found that the retrieved extinction coefficients at 17 nm at heights between 280 and 450 km were larger by up to 30–60% in comparison with those in the quiet geomagnetic state. Below 280 km, the extinction remained the same during the storm. The DTM-2013 model describes this variation of extinction during the storm very well, whereas the NRLMSISE-00 model failed to reproduce it.

In conclusion, our investigation has shown that the method of solar occultations in the EUV range is a very sensitive tool to study the properties of the upper atmosphere at heights 200–500 km and its response to solar and geomagnetic activity.

Acknowledgments

The research leading to these results has received funding from the European Commission’s Seventh Framework Programme (FP7/2007-2013) under the Grant Agreement eHeroes (project No 284461, http://www.eheroes.eu). SWAP is a project of the Centre Spatial de Liège and the Royal Observatory of Belgium funded by the Belgian Federal Science Policy Office (BELSPO). The SPIRIT and TESIS experiments have been developed and carried out by the Laboratory of X-ray Astronomy of the Sun, P.N. Lebedev Physical Institute of the Russian Academy of Sciences, under the CORONAS International solar space program. The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

Annex A: Retrieval of extinction coefficients from measured profiles

A1. Reduction of Sun measurements to a point source

Let us consider the occultation scheme (Fig. 1) and rewrite Eq. (1) with notations more convenient for the analysis. The received solar radiant intensity I^m (r₀, y₀) depends on the ray perigee position of the center of the solar disk r₀ and it does not depend on the transversal position (to the x-r plane) y₀ of the Sun. This intensity is an integral over the disk surface D of the surface density of the radiant intensity J(r, y) (the radiant intensity of the received emission per unit disk surface). The value of J is also independent of the transversal position y of the corresponding element of the solar surface. Taking this into account, one has$$I^{\mathrm{m}}\left(r_0,y_0\right)=\underset{D}{\overset{}{\iint }}J\left(r,y\right)\mathrm{d}r\mathrm{d}y=\underset{r_0-r_{\mathrm{s}}}{\overset{r_0+r_{\mathrm{s}}}{\int }}J\left(r\right)\mathrm{d}r\underset{y_0-\sqrt{r_{\mathrm{s}}^2-{(r-r_0)}^2}}{\overset{y_0+\sqrt{r_{\mathrm{s}}^2-{(r-r_0)}^2}}{\int }}\mathrm{d}y,$$(A1)where r_s is the radius of the Sun projection around the tangential point (i.e., r_s = dh in Fig. 1).

This expression is valid for a uniformly radiating solar disk. Multiplying J(r) by the disk area, one obtains the intensity dependence I(r), such as it would be for the whole Sun located at the ray perigee position r. The relationship between the measured dependence I^m(r₀) and I(r) is expressed by the integral equation:$$I^{\mathrm{m}}\left(r_0\right)=\underset{r_0-r_{\mathrm{s}}}{\overset{r_0+r_{\mathrm{s}}}{\int }}I\left(r\right)\frac{2}{\pi r_{\mathrm{s}}^2}\sqrt{r_{\mathrm{s}}^2-{(r-r_0)}^2}\mathrm{d}r=\underset{-\infty }{\overset{\infty }{\int }}I\left(r\right)F(r-r_0)\mathrm{d}r,$$(A2) $$F=\left\{\begin{array}{cc}\frac{2}{\pi r_{\mathrm{s}}^2}\sqrt{r_{\mathrm{s}}^2-{(r-r_0)}^2},& |r-r_0|\le r_{\mathrm{s}}\\ 0,& |r-r_0|>r_{\mathrm{s}}\end{array}\right\}.$$

The kernel of (A2) is normalized to unity, hence in the case of slow variation of I(r) over the solar disk one has I^m(r₀) ≈ I(r). The solution of Eq. (A2) gives the dependence of radiance intensity I(r) of an equivalent point source located at a fixed ray perigee r. It is a Volterra-kind equation with variable integration limits, and it is known as a weakly ill-posed problem. This equation can also be expressed as a Fredholm integral equation of the 1-st kind with the kernel F (Eq. (A2)) that can easily be solved using Tikhonov’s method of generalized discrepancy. In the case of an inhomogeneous distribution of solar radiation over the disk, this distribution should be taken into account in the formation of the kernel F. The retrieved dependence I(r) satisfies the evident relationship following from (1):$$I\left(r\right)=I_0\exp \left(-\tau \left(r\right)\right),$$(A3)where τ is the optical depth along the ray of perigee r and I₀ is the solar radiant intensity outside the atmosphere. From (A3) one can express the difference of the optical depths corresponding to the rays with perigee r and r_m as$$\tau \left(r_{\mathrm{m}}\right)-\tau \left(r\right)=\ln \left(\frac{I\left(r\right)}{I\left(r_{\mathrm{m}}\right)}\right).$$(A4)

Under the assumption that the Earth atmosphere is spherically symmetric, this difference is determined by the absorption altitude profile γ(r) as:$$\Delta \tau \left(r_0\right)=\tau \left(r_{\mathrm{m}}\right)-\tau \left(r_0\right)=\underset{r_0}{\overset{\infty }{\int }}\gamma \left(r\right)\Delta K(r,r_0)\mathrm{d}r=\underset{r_0}{\overset{\infty }{\int }}\gamma \left(r\right)\left[K(r,r_{\mathrm{m}})-K(r,r_0)\right]\mathrm{d}r,$$(A5) $$K\left(r,r_0\right)=\left\{\begin{array}{cc}\frac{2r}{\sqrt{r^2-r_0^2}},& r<r_H\\ \frac{r}{\sqrt{r^2-r_0^2}},& r\ge r_H\end{array}\right\},$$where r_H is the radial distance of the receiver.

In our case, we can restrict the analysis to the situation where r < r_H. In this case, the corresponding inverse problems can be solved using the Abel inverse transform (in the refraction inverse problem) or the inversion of the Volterra integral equation of the 1-st kind in its classical statement (Gaikovich 2004). If there would be some atmospheric contribution above the level r = r_H, Eq. (A5) transforms into a Volterra integral equation of the 1-st kind in the ill-posed statement (Gaikovich 2004, p. 21).

A high-accuracy solution (like that demonstrated in this paper) can be achieved for the layers below the height of the receiver since the contribution of the perigee vicinity is strongly dominant because of the weak peculiarity in the kernel of (A5).

All our results were obtained by data at negative elevation angles (as illustrated in Fig. 1). It is nevertheless worth noting that by performing measurements also at the respective symmetric positive angles and looking at the difference between the two, one can guarantee that the solution does not depend on the atmosphere above the height of the receiver. In this case, instead of (A5), one has an ordinary Volterra equation of the 1-st kind with a weak peculiarity that can easily be solved (Gaikovich 2004):$$\Delta \tau \left(r_0\right)={\tau }^{-}\left(r_0\right)-{\tau }^{+}\left(r_0\right)=\underset{r_0}{\overset{r_Y}{\int }}\gamma \left(r\right)\Delta K(r,r_0)\mathrm{d}r=\underset{r_0}{\overset{r_Y}{\int }}\gamma \left(r\right)\left[K^{-}\left(r,r_0\right)-K^{+}(r,r_0)\right]\mathrm{d}r,$$(A6) $$\Delta K\left(r,r_0\right)=\frac{2r}{\sqrt{r^2-r_0^2}}.$$

In general, we have to deal with the much more difficult problem of finding the solution of (A5) both below and above the altitude of the receiver. A previous application of the described method to the reconstruction of the refraction parameters of the atmosphere shows that at altitudes much lower than the receiver altitude one has a sufficiently accurate solution (Gaikovich et al. 1983; Gaikovich 2004); but at heights comparable to and higher than the receiver altitude, the kernel of (A5) changes – it has no weak peculiarity above this level. Hence, the solution quality in this range is inherent to ill-posed problems related to Fredholm integral equations of the 1-st kind: the solution has a weak convergence when reducing the measurement errors, and the retrieval accuracy is sharply reduced with the increase of the altitude. An estimation of the retrieval accuracy can only be obtained, as in any ill-posed problem, from numerical simulation taking into consideration the convergence properties of the regularized solution.

A2. Method of absorption profiling and results of the numerical simulation

The considered problem consists in solving the integral Eqs. (A2) and (A5). Both equations can also be expressed as an integral Fredholm equation of the 1-st kind$$z_{\mathrm{m}}\left(r_0\right)=\underset{r_1}{\overset{r_2}{\int }}y\left(r\right)K(r-r_0)\mathrm{d}r=K_{h}y,$$(A7)where K_h is the approximation of the equation kernel K. The upper limit should be chosen so that the contribution of the range above this limit is smaller than the measurement errors. The lower limit can be chosen arbitrarily in the range r_0min – r_s ≤ r₁ < r₂, where r_0min is the lower value of the perigee of the center of the solar disk in measurements. This makes it possible to solve the problem with various values of r1. This property can be useful when measurement errors are strongly altitude-dependent.

To solve Eq. (A7), it is necessary to use regularization methods which involve additional a priori information about the function to be retrieved. The Tikhonov’s method of generalized discrepancy considered here (Tikhonov & Arsenin 1977) uses the information that the solution belongs to the set of square-integrable functions with square-integrable derivatives. In this method, an approximate solution y^α of (A3) is found using the parameter of measurement errors δ and the parameter of the kernel error h_k that satisfy the conditions:$${\left|\right|z-z_{\mathrm{m}}\left|\right|}_{L_2}^2\le {\delta }^2, \hspace{1em}{\left|\right|K-K_h\left|\right|}_{W_1^2\to L_2}\le h_{\mathrm{k}},$$(A8)where z is the left-hand side of (A7) that corresponds to the exact solution y; z_m represents the measurements (affected by measurement errors). According to Tikhonov’s method, the approximate solution of (A7) is found from minimizing the generalized discrepancy:$$M_{\alpha }\left[y^{\alpha }\right]={\left|\right|K_h{y}^{\alpha }-z_{\mathrm{m}}\left|\right|}_{L_2}^2+\alpha {\left|\right|y^{\alpha }\left|\right|}_{W_2^2}^2$$(A9)with the additional condition (principle of generalized discrepancy)$${\left|\right|K_h{y}^{\alpha }-z_{\mathrm{m}}\left|\right|}_{L_2}^2={\left(\delta +h_{\mathrm{k}}{\left|\right|y^{\alpha }\left|\right|}_{w_1^2}\right)}^2,$$(A10)and makes it possible to find the value of the regularization parameter at given values of δ and h.

An important advantage of this method is the convergence of the approximate solution to the exact one in the metric W₁², and, hence, the uniform convergence in the metric C (with the maximum modulus as a norm) at (δ, h_k) → 0 estimated in the integral metric L₂. So there is no need to use large statistical ensembles to obtain representative estimations of the retrieval accuracy as it is necessary in other methods with integral or mean square convergence. It is sufficient to run the numerical simulation for typical and, maybe, for extreme exact solutions to estimate the maximum error of the retrieved distribution. There is some uncertainty in the measurement error estimation related to the random origin of errors that should be taken into account. It is of advantage to search for a solution as the deviation from a reasonable model function. This deviation is typically a smoother and simpler function, and such functions are much better retrieved in solving ill-posed problems.

The accuracy of the proposed method has been studied by numerical simulations. These simulations involved the following steps: (a) for an initial simulated distribution γ(r) the exact left-hand side of (A3) I(r) is calculated; (b) random Gauss-distributed errors with rms that correspond to the real measurement accuracy are added to I(r) to obtain “measurement data”; (c) the inverse problem (A3) is solved to obtain Δτ(r₀); (d) these data are used to retrieve the distribution γ(r) from the solution of the inverse problem (A6); (e) the retrieved distribution γ(r) is compared to the initial one.

In numerical simulations, we used the initial distributions γ(h) obtained from deviations from the atmosphere model:$$\gamma \left(h\right)={\gamma }_0\left(h\right)+A\exp \left[-{\left(\frac{h-h_{\max }}{\Delta h}\right)}^2\right].$$(A11)

Where the atmospheric model$${\gamma }_0\left(h\right)={10}^{(k_1+k_{2}h+k_3{h}^2)},$$(A12)is a function of h the altitude above the earth surface, as well as of several coefficients here obtained from the NRLMSISE-00 model: k₁ = −5.94, k₂ = −0.00797, k₃ = 2.36 × 10⁻⁷ for wavelength 17.5 nm.

The measure of retrieval error in solving (A6) $${\delta }^2=\sup {\left|\left|\mathbf{K}\gamma -{\tau }^{\delta }\right|\right|}_{L_2}^2=\frac{1}{\Delta r}\underset{\Delta r}{\overset{}{\int }}{\left[\tau \left(r\right)-{\tau }^{\delta }\left(r\right)\right]}^2\mathrm{d}r,$$(A13)was estimated from measurements as δ = 0.05.

From results of the above-described numerical simulation with the level of random errors δ = 0.05 (in the integral metric L₂), true estimations of the accuracy of the retrieved profiles γ(h) were obtained in the maximum modulus metric. They demonstrate that the accuracy of retrieval depends on the depth scale ∆h of the profile perturbation: it changes from 10% for large-scale variations at ∆h ~ 50 km or more up to 30% for small-scale variations at ∆h ~ 10 km.

All Tables

All Figures

	Figure 1. Scheme of the occultation measurements. R0 – radius of the Earth, Horb – altitude of the satellite orbit. Other notations are given in the text.
In the text

	Figure 2. (A) Sunspot numbers in solar cycles 23 and 24 and observation times of the SPIRIT, TESIS, and SWAP data (explanations in the text); (B) transmittance profiles for SPIRIT (solar maximum, data set 1, red line), TESIS (solar minimum, black line), and SWAP (rising phase of cycle 24, blue line).
In the text

	Figure 3. Map of occultation locations: red symbols – locations of sub-satellite points for SPIRIT (data set 1), TESIS and SWAP, blue symbols – locations of sub-LOS points.
In the text

	Figure 4. Height profiles of extinction coefficients of the main atmospheric constituents calculated with the NRLMSISE-00 and DTM-2013 models for the EUV wavelength 17 nm in the quiet geomagnetic state. (A) At solar maximum (the SPIRIT data set 1), (B) at solar minimum (the TESIS data set).
In the text

Figure 5. Comparison of the transmittance height profiles, presented in Figure 2B, with those calculated by the NRLMSISE-00 and DTM-2013 models. The black lines display the averaged measurement points (for TESIS individual points are also shown), the red/blue lines show simulations by the NRLMSISE-00 and DTM-2013 models. In the SPIRIT and SWAP cases, the error bars correspond to the transmittance determination errors, in the TESIS case to the dispersion of data averaged over multiple profiles.

In the text

	Figure 6. Comparison of the modeled and retrieved extinction coefficients as a function of height for the SPIRIT data set 1 (March 5, 2002, left) and the SWAP data set (December 22, 2010, right). Black line – the data retrieved from measurements, red/blue lines – calculations by the NRLMSISE-00 and DTM-2013 models.
In the text

	Figure 7. Variation of the retrieved and modeled extinction coefficients during the development of the geomagnetic storm of 23 May 2002: (A) before the storm (the SPIRIT data set 2, Dst = −1 nT), (B) 15 h later (the SPIRIT data set 3, Dst = −75 nT). The line colors are the same as in Figure 6.
In the text