An inverse model to correct for the effects of post-depositional processing on ice-core nitrate and its isotopes: model framework and applications at Summit, Greenland, and Dome C, Antarctica

Comprehensive evaluation of the effects of post-depositional processing is a prerequisite for appropriately interpreting ice-core records of nitrate concentration and isotopes. In this study, we developed an inverse model that uses archived snow/ice-core nitrate signals to reconstruct primary nitrate flux (i.e., the deposition flux of nitrate to surface snow that originates from long-range transport or stratospheric input) and its isotopes (δ¹⁵N and Δ¹⁷O). The model was then applied to two polar sites, Summit, Greenland, and Dome C, Antarctica, using measured snowpack nitrate concentration and isotope profiles in the top few meters. At Summit, the model successfully reproduced the observed atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) and their seasonality. The model was also able to reasonably reproduce the observed snowpack nitrate profiles at Dome C as well as the skin layer and atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) at the annual scale. The calculated F_pri at Summit was 6.9 × 10⁻⁶ kgN m² a⁻¹, and the calculated Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri is consistent with atmospheric observations in the Northern Hemisphere. However, the calculated δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) of F_pri displays an opposite seasonal pattern to atmospheric observations in the northern mid-latitudes, but it is consistent with observations in two Arctic coastal sites. The calculated F_pri at Dome C varies from 1.5 to 2.2 × 10⁻⁶ kgN m⁻² a⁻¹, with δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) of F_pri varying from 6.2 ‰ to 29.3 ‰ and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri varying from 48.8 ‰ to 52.6 ‰. The calculated F_pri at Dome C is close to the previous estimated stratospheric denitrification flux in Antarctica, and the high δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri at Dome C also point towards the dominant role of stratospheric origin of primary nitrate to Dome C.

1 Introduction

Nitrate ion (NO${}_{\mathrm{3}}^{-}$) is routinely measured in polar snow and ice cores. The precursor of atmospheric nitrate is nitrogen oxide NO_x (= NO + NO₂), which plays a fundamental role in the production of tropospheric ozone and interconversion of atmospheric HO_x (= OH + HO₂) radicals (Seinfeld and Pandis, 2016; Sillman, 1999). Given the potential link between ice-core nitrate and atmospheric NO_x, some previous studies proposed that ice-core nitrate records could be used to derive information regarding past atmospheric NO_x abundance (Dibb et al., 1998; Röthlisberger et al., 2000). In addition, the oxygen isotope mass-independent fractionation signal (Δ¹⁷O = δ¹⁷O − 0.52 × δ¹⁸O) of nitrate is a reliable proxy of atmospheric ${\mathrm{O}}_{\mathrm{3}}/{\mathrm{HO}}_x$ ratio and is directly related to atmospheric oxidizing environment (Alexander et al., 2004; Alexander and Mickley, 2015; Geng et al., 2017; Sofen et al., 2014). These unique features render ice-core nitrate a potentially useful proxy to retrieve information on atmospheric oxidation environment in the past (Alexander and Mickley, 2015).

Interpretations of ice-core nitrate records are, however, not straightforward (Wolff et al., 2008). Unlike other less reactive species in ice cores such as sulfate, ice-core nitrate may not be able to directly track its atmospheric abundance (Iizuka et al., 2018). To link ice-core nitrate to atmospheric NO_x abundance, other information including the conversion rate of NO_x to nitrate, the mean lifetime of atmospheric nitrate, and the impact of post-depositional processing must be considered (Wolff, 1995; Wolff et al., 2008). Among these factors, the post-depositional processing of snow nitrate is the first gap in linking ice-core nitrate to atmospheric nitrate and/or NO_x.

Snow nitrate is reactive under exposure to sunlight and can be photolyzed to form NO_x and HONO (Honrath et al., 2002; Chu and Anastasio, 2003), which is rapidly transported to the overlying atmosphere via diffusion and convection (Zatko et al., 2013). These photoproducts subsequently reform nitrate (we further refer to this as snow-sourced nitrate) which is redeposited locally or exported away, leading to a recycling of nitrate at the air–snow interface (Erbland et al., 2013; Frey et al., 2009). The reformed nitrate would inherit Δ¹⁷O signals under local oxidation conditions that is different from primary nitrate, and the re-deposition of atmospheric nitrate could also result in nitrogen isotopic fractionation depending on the different deposition mechanisms (Erbland et al., 2013; Jiang et al., 2022). This post-depositional processing not only disturbs the link between nitrate in snow and its atmospheric precursors but also alters its isotopic signals (Erbland et al., 2013; Jiang et al., 2021, 2022; Shi et al., 2015). But since these processes are initiated by sunlight, the post-depositional processing is muted in polar winter when sunlight is absent.

It is expected that the degree of post-depositional processing varies with changes in factors such as snow accumulation rate under different climates (Akers et al., 2022b; Geng et al., 2015), causing corresponding shifts in the preserved nitrate signals. For example, the lower snow accumulation rate in glacial times would favor a higher degree of post-depositional processing with elevated δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) relative to the Holocene as reflected by the GISP2 ice-core records (Geng et al., 2015; Hastings et al., 2005). Moreover, both observational and modeling studies have suggested that at sites with relatively high snow accumulation rates such as Summit, Greenland, the post-depositional processing of snow nitrate under present-day conditions also has a significant impact on seasonal δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) variations, although its integral effects at the annual scale are limited (Jiang et al., 2021, 2022). In addition, the Δ¹⁷O of snow nitrate would also be altered via secondary chemistry during photolysis on snow grain (i.e., the cage effect), and this effect is enhanced with lower snow accumulation rates (Erbland et al., 2013; Frey et al., 2009; McCabe et al., 2005; Meusinger et al., 2014). Thus, it is critical to evaluate the impact of post-depositional processing on ice-core nitrate records before interpretation, especially for records covering different climates with changes in snow accumulation rates.

Primary nitrate to the polar ice sheets mainly originates from midlatitudes via long-range transport and with extra contributions from stratospheric input (Lee et al., 2014; Legrand and Delmas, 1986; Fischer et al., 1998; Savarino et al., 2007). To build the link between ice-core and atmospheric nitrate, Geng et al. (2015) proposed a simple method of using δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) to estimate the fractional loss of snow nitrate caused by post-depositional processing. This method takes advantage of the high sensitivity of δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) to the degree of photolytic loss (Erbland et al., 2013; Frey et al., 2009). If δ¹⁵N of the initially deposited nitrate can be assumed, the residual fraction of snow nitrate can be calculated by applying a Rayleigh type isotope fractionation model. The photolysis fractionation constant (¹⁵ε_p) can be estimated via the prescribed actinic flux spectrum and the absorption cross section for different nitrate isotopologues (Berhanu et al., 2014). Based on this method, Geng et al. (2015) estimated that as much as 45 %–53 % of snow nitrate was lost after deposition during the last glacial time in the GISP2 ice-core record. However, it is difficult to justify the assumed δ¹⁵N of deposited nitrate under different climates, and the method cannot correct for post-depositional modification of Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$).

Erbland et al. (2015) developed a 1-D snow photochemistry model (TRANSITS, https://github.com/JZxxhh/TRANSITS-model, last access: 29 January 2024) that quantifies the effects of post-depositional processing on the preservations of nitrate and its isotopes in ice cores. The model comprises a series of physicochemical processes, including UV photolysis of snow nitrate, emission of NO_x to the overlying atmosphere, local oxidation, and nitrate deposition. In addition, changes in the isotopic composition of nitrate (δ¹⁵N and Δ¹⁷O) at each step of the post-depositional processing are also explicitly incorporated. Recently, Shi et al. (2023) extended or followed the TRANSITS model framework to include snowpack δ¹⁸O(NO${}_{\mathrm{3}}^{-}$) simulation during the preservation of nitrate in snow. The latter was built upon the same chemical processes related to modeling Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) changes during the post-depositional processing. However, the fractionation factor of δ¹⁸O during snow nitrate photolysis (¹⁸ε_p) had to be scaled to reproduce the observations. In this case it remains unclear why the theoretical fractionation factor calculated using the ZPE (zero-point energy) shifted method (Frey et al., 2009) works well on δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) but not on δ¹⁸O(NO${}_{\mathrm{3}}^{-}$). Nevertheless, the uncertainties associated with δ¹⁸O fractionations during snow nitrate photolysis and other processes (e.g., the cage effect, reformation of nitrate from NO₂) make this simulation less useful and reliable than for Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$), for which there are much fewer influencing factors and are easier to constrain. The TRANSITS model has been applied in various locations with different snow accumulation rates and well reproduced the observed snowpack nitrate and isotope profiles (Erbland et al., 2015; Jiang et al., 2021; Winton et al., 2020; Zatko et al., 2016). Based on model sensitivity tests, Erbland et al. (2015) proposed a framework to correct for the effects of post-depositional processing and to retrieve atmospheric information related to F_pri at Dome C. However, the framework is rather complicated, and it assumes δ¹⁵N of the archived nitrate is exclusively determined by the degree of nitrate post-depositional processing. Therefore, the framework cannot be applied to sites with moderate or high snow accumulation rates such as WAIS Divide, Antarctica, and Summit, Greenland, where factors other than post-depositional processing may also contribute to δ¹⁵N variations across different periods and/or climates (Hastings et al., 2005; Jiang et al., 2021).

In summary, TRANSITS is a forward model, and it requires prior knowledge of the distribution (e.g., weekly or monthly) of primary nitrate flux and isotopes as model inputs, which is usually unavailable due to the lack of direct observations. In this study, we developed an inverse modeling framework (i.e., the inverse of the TRANSITS model) that uses snowpack and/or ice-core preserved nitrate signals (concentrations and isotopes) as model inputs, and properties of primary nitrate including its flux and isotopes (δ¹⁵N and Δ¹⁷O) can be directly retrieved with constraints from snow accumulation rate and other known parameters (e.g., snow physicochemical properties). We assessed the model with observations at Summit, Greenland, and Dome C, Antarctica, two representative sites with approximately the high-end and low-end snow accumulation rates at present-day conditions.

2 Model description

The inverse model is designed based on the framework of the TRANSITS model but in an opposite direction of operating flows. The principle of the inverse model is that the archived snow nitrate concentration and isotope profiles from measurements are treated as model input, and they evolve inversely over time through the snow photic zone (defined as 3 times of the snow e-folding depth where the radiation decreases to $\mathrm{1}/e$ of its initial intensity at snow surface) to recover their initial states at the time of deposition, thus providing the initial isotope compositions and deposition fluxes before being affected by any post-depositional effects. The primary nitrate flux and its isotopes can be further obtained by solving the mass balance equations in the atmosphere box. A schematic view of the inverse model is shown in Fig. 1 with arrows pointing toward the model direction flow (i.e. inverse of the real physical processes). Major parameters in the inverse model and their definitions are listed in Table 1.

Table 1 List of major parameters used in the inverse model.

^* Three types of light absorption impurity are considered in the inverse mode: black carbon, mineral dust, and organic humic-like substance (HULIS).

Figure 1 Schematic plot of the model domains of the inverse model including the atmospheric box, the snow photic zone and the archived snow layers, where f_exp represents the fraction of nitrate exported from the site of photolysis. The nitrate isotopic and mass balance relationships on snow grains during photolysis are also shown, where f_p represents the fraction of snow nitrate being photolyzed, and f_c represents the fraction of photolyzed nitrate experiencing the cage effect (i.e., exchange of oxygen isotopes with snow water). FA represents the archived nitrate flux.

The inverse model inherits most of the original processes and features in TRANSITS but with several modifications. In accordance with the TRANSITS model, the domains of the inverse model are represented by a 1-D atmosphere and snow column. As shown in Fig. 1, the model contains three vertical parts, including the overlying atmospheric boundary layer, which is treated as a single well-mixed box, the underlying snowpack, which is further separated into a snow photic zone, and the archived snow layers beneath the photic zone. The model time step is set to be 1 week by default. During each time step, the mass conservation equations in the atmospheric box are represented as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{m}_{\mathrm{a}}}{\mathrm{d}t}}=F_{\mathrm{pri}}+\mathrm{FP}-\mathrm{FE}-\mathrm{FD}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}(m_{\mathrm{a}}\times {\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{a}})}{\mathrm{d}t}}& =F_{\mathrm{pri}}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}\left(F_{\mathrm{pri}}\right)+\mathrm{FP}\times \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FP}\right)\\ & -\mathrm{FE}\times \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FP}\right)-\mathrm{FD}\times \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FD}\right)\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}(m_{\mathrm{a}}\times \mathrm{\Delta }{}^{\mathrm{17}}{\mathrm{O}}_{\mathrm{a}})}{\mathrm{d}t}}& =F_{\mathrm{pri}}\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(F_{\mathrm{pri}}\right)+\mathrm{FP}\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FP}\right)\\ & -\mathrm{FE}\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FE}\right)-\mathrm{FD}\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FD}\right),\end{array}\end{array}$$

where the subscript “a” represents the atmospheric box (i.e., m_a refers to the mass of atmospheric nitrate), and δ¹⁵N_a and Δ¹⁷O_a refer to δ¹⁵N and Δ¹⁷O of atmospheric nitrate, respectively. Different nitrate fluxes transported in and out of the atmospheric box are denoted as FP, FE, and FD, where FP refers to the photolytic nitrate flux (the snow-sourced nitrate), FD refers to the atmospheric deposition nitrate flux, and FE refers to the exported nitrate flux that is horizontally transported out of the atmospheric box via air flow. Following Erbland et al. (2015), FE is assumed to be a portion (f_exp) of FP (i.e., FE = f_exp×FP) and maintains the isotopic signatures of FP.

In Eqs. (1)–(3), the LHS (left-hand side) terms are 2–3 orders of magnitude smaller than nitrate fluxes in and out of the atmospheric box. Erbland et al. (2015) showed that the atmospheric nitrate mass was a factor of ∼ 10⁻³ smaller than the surface snow nitrate reservoir at Dome C, and similar results were also found at Summit in Jiang et al. (2021). Thus, $\mathrm{d}\left(x\right)/\mathrm{d}t$ is assumed to be zero at each time step (i.e., species and isotope compositions in the atmosphere are considered at steady state), which leads to simplified formulas for calculating F_pri via Eqs. (4)–(6) as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{F}_{\mathrm{pri}}\approx \mathrm{FD}-\mathrm{FP}\left(\mathrm{1}-f_{\exp }\right)\text{(5)}& {\displaystyle }\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(F_{\mathrm{pri}}\right)\approx {\displaystyle \frac{\begin{array}{c}\mathrm{FD}\times \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FD}\right)\\ -\mathrm{FP}\left(\mathrm{1}-f_{\exp }\right)\times \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FP}\right)\end{array}}{\mathrm{FD}-\mathrm{FP}\left(\mathrm{1}-f_{\exp }\right)}}\text{(6)}& {\displaystyle }\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(F_{\mathrm{pri}}\right)\approx {\displaystyle \frac{\begin{array}{c}\mathrm{FD}\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FD}\right)\\ -\mathrm{FP}\left(\mathrm{1}-f_{\exp }\right)\times \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FP}\right)\end{array}}{\mathrm{FD}-\mathrm{FP}\left(\mathrm{1}-f_{\exp }\right)}}.\end{array}$$

Hence, if the magnitude and isotopic compositions of FP and FD in each time step are known, F_pri can be calculated. FP and FD calculated from the inverse evolution of snowpack nitrate are described in the following sections.

2.1 The backward evolution of snowpack nitrate

Starting with an arbitrary snowpack nitrate depth profile at a given time step, changes in nitrate concentration and isotopic compositions (δ¹⁵N and Δ¹⁷O) in a certain snow layer in the photic zone induced by photolysis can be calculated as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }c\left({\mathrm{SN}}_n^{\prime }\right)={\displaystyle \frac{c\left({\mathrm{SN}}_n\right)}{(\mathrm{1}-f_{\mathrm{p}})+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\text{(8)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_n^{\prime }\right)& =\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_n\right)\\ & -{\displaystyle \frac{(\mathrm{1}-f_{\mathrm{p}})(\mathrm{1}-f_{\mathrm{c}})\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}\ln (\mathrm{1}-f_{\mathrm{p}})}{(\mathrm{1}-f_{\mathrm{p}})+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\end{array}\text{(9)}& {\displaystyle }\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}_n^{\prime }\right)=\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}_n\right){\displaystyle \frac{(\mathrm{1}-f_{\mathrm{p}})+f_{\mathrm{c}}{f}_{\mathrm{p}}}{(\mathrm{1}-f_{\mathrm{p}})+\frac{\mathrm{2}}{\mathrm{3}}{f}_{\mathrm{c}}{f}_{\mathrm{p}}}},\end{array}$$

where c represents the nitrate concentration and SN_n refers to the nth snowpack layer, respectively, and the quotation mark in superscript refers to the initial state before being photolyzed at each time step. These equations are based on the nitrate mass and isotopic balances on snow grains during photolysis as shown in Fig. 1, and detailed derivations of these equations can be found in Appendix A.

In Eqs. (7)–(9), f_p represents the fraction of snow nitrate that undergoes photolysis at each time step, and f_c represents the fraction of nitrate photolysis intermediate undergoing the cage effect (Meusinger et al., 2014), which leads to apparent oxygen isotope exchange with water and lowers Δ¹⁷O by a factor of $\mathrm{2}/\mathrm{3}$. The potential isotope effect on δ¹⁵N during cage effect remains unknown and is not considered. The value of f_p is calculated by the first-order reaction of nitrate photolysis:

$$\begin{array}{}\text{(10)}& f_{\mathrm{p}}=\mathrm{1}-\exp \left(-\underset{\mathrm{0}}{\overset{\mathrm{d}t}{\int }}J(t,z)\mathrm{d}t\right),\end{array}$$

where J represents the rate constant of nitrate photolysis that varies with time and depth of the snow layer. J is calculated from actinic flux (I), the quantum yield (Φ), and the absorption cross section (σ) of nitrate photolysis as follows:

$$\begin{array}{}\text{(11)}& J\left(t,z\right)=\underset{\mathrm{280}\mathrm{nm}}{\overset{\mathrm{350}\mathrm{nm}}{\int }}\mathrm{\Phi }\left(\mathit{\lambda }\right)\times {\mathit{\sigma }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\left(\mathit{\lambda }\right)\times I\left(z,\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }.\end{array}$$

The rate constant of ¹⁵NO${}_{\mathrm{3}}^{-}$ photolysis (J^∗) is also calculated from the absorption cross section of the heavy isotopologue from Berhanu et al. (2014), and the photolysis fractionation constant for nitrogen isotope ε_p is calculated via

$$\begin{array}{}\text{(12)}& {\mathit{\epsilon }}_{\mathrm{p}}\left(t,z\right)={\displaystyle \frac{J^{\ast }\left(t,z\right)}{J\left(t,z\right)}}-\mathrm{1}.\end{array}$$

The solar zenith angle changes with time during each time step, leading to changes in the spectrum of actinic flux and subsequently changes in ε_p. To simplify the calculation, in Eq. (8) $\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}$ in a certain week is calculated by the weighted average of nitrogen isotope fractionation constant over the durations of different solar zenith angles (0–90°). The radiative transfer in snowpack is calculated using the parameterization from Zatko et al. (2013) to achieve fast online calculations, and this parameterization has been shown to be capable of providing consistent results with a high-order snowpack radiative transfer model DISORT (Zatko et al., 2013). The upper boundary conditions for the parameterization, i.e., the direct and diffuse components of the irradiance at the snow surface, are calculated offline using the Troposphere Ultraviolet and Visible (TUV) radiation model (Madronich et al., 1998) at different total column ozone (TCO) and solar zenith angle conditions.

The relationships between c(SN_n), δ¹⁵N(SN_n), and Δ¹⁷O(SN_n) and c(SN${}_n^{\prime }$), δ¹⁵N(SN${}_n^{\prime }$), and Δ¹⁷O(SN${}_n^{\prime }$) in the snowpack are illustrated in Fig. 2. Here, c(SN_n), δ¹⁵N(SN_n), and Δ¹⁷O(SN_n) represent the values after photolysis in the nth layer at a certain time step, while c(SN${}_n^{\prime }$), δ¹⁵N(SN${}_n^{\prime }$), and Δ¹⁷O(SN${}_n^{\prime }$) denote the values before photolysis at this time step, which are calculated using Eqs. 7–9. c(SN${}_n^{\prime }$), δ¹⁵N(SN${}_n^{\prime }$), and Δ¹⁷O(SN${}_n^{\prime }$) also correspond to the values after photolysis in the prior time step when it was in the (n−1)th layer. By repeating this operation, the initially deposited values of nitrate concentration and isotopes for a given snow layer without influence from the photo-driven post-depositional processing (i.e., when this layer was at the surface) can be calculated, which is be further linked to FD.

Figure 2 Schematic of the evolution of snowpack nitrate from the archived layer to its initial state at the snow surface. The quotation mark in the superscript of the bracket represents the status of snow nitrate before photolysis at each time step.

2.2 Determinations of FP and FD

FP and FD are determined during the inverse evolution of snowpack nitrate profiles. As shown in Fig. 1, photolytic nitrate flux and its δ¹⁵N from the nth snow layer can be calculated via the mass balance relationships:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{FP}}_n=c\left({\mathrm{SN}}_n^{\prime }\right)f_{\mathrm{p}}\left(\mathrm{1}-f_{\mathrm{c}}\right){\mathit{\rho }}_{\mathrm{snow}}{d}_n/\mathrm{\Delta }t\text{(14)}& {\displaystyle }\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{FP}}_n\right)=\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_n\right)-{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}(\mathrm{1}-f_{\mathrm{p}})\ln (\mathrm{1}-f_{\mathrm{p}})}{f_{\mathrm{p}}}}.\end{array}$$

In Eq. (13), ρ_snow is the density of snow, d_n is the thickness of the nth snow layer, which is equal to the accumulated snow thickness at one time step, and Δt is the default model time step (1 week). Equation (14) implicitly assumes that the reformed nitrate in the overlying atmosphere keeps the same δ¹⁵N signals of the snow-emitted photoproduct of NO₂ because of isotope mass balance; i.e., essentially all NO₂ is oxidized into nitrate at one time step. FP emitted from the whole snowpack and its δ¹⁵N can be calculated by

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{FP}=\sum {\mathrm{FP}}_n\text{(16)}& {\displaystyle }\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FP}\right)={\displaystyle \frac{\sum {\mathrm{FP}}_n\cdot \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{FP}}_n\right)}{\sum {\mathrm{FP}}_n}}.\end{array}$$

For Δ¹⁷O of FP, extra knowledge of the oxidizing agent concentrations in the local atmosphere including HO₂, RO₂, and O₃ must be provided (Appendix B). This is because the emitted NO_x would achieve photochemical steady state rapidly, thus erasing any original Δ¹⁷O signal inherited from the snowpack nitrate. During the subsequent oxidation of atmospheric NO₂, one more oxygen atom inherited from the oxidants (e.g., OH or BrO) is incorporated into one newly formed HNO₃ molecule. Thus, Δ¹⁷(FP) can be represented by $\mathrm{2}/\mathrm{3}$ of Δ¹⁷O(NO₂) plus $\mathrm{1}/\mathrm{3}$ of Δ¹⁷O(oxidant).

FD and its isotopic signals can be obtained from the uppermost snow layer before photolysis occurs as illustrated in Fig. 2:

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathrm{FD}=c\left({\mathrm{SN}}_{\mathrm{0}}^{\prime }\right){\mathit{\rho }}_{\mathrm{snow}}{d}_{\mathrm{0}}/\mathrm{\Delta }t\text{(18)}& {\displaystyle }\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{FD}\right)=\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_{\mathrm{0}}^{\prime }\right)\text{(19)}& {\displaystyle }\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{FD}\right)=\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}_{\mathrm{0}}^{\prime }\right).\end{array}$$

The calculated FP and FD in each time step are further used to calculate F_pri according to Eqs. (4)–(6).

2.3 The choice of model initial conditions

To run the model, an appropriate archival snow nitrate profile with known concentration and isotopic composition (δ¹⁵N and Δ¹⁷O) should be assigned as model initial conditions with seasonal or monthly resolution, though ideally weekly or finer-resolution data are the best. The archived nitrate profile could be dated by using various types of seasonal markers, such as the δ¹⁸O of H₂O, the ion concentrations or their ratios, and the snow accumulation rates (Hastings et al., 2004; Furukawa et al., 2017; Dibb et al., 2007). As long as the archived snow nitrate profiles (i.e., snow nitrate concentration and isotopes below the photic zone) are given, the model can calculate nitrate concentrations and isotopes throughout the photic zone and those in the atmosphere. The latter is considered as the atmospheric signals before being affected by post-depositional processing.

3 Model evaluations

Because we lack direct observations of primary nitrate, we evaluated the model performance with other kinds of observations, including nitrate isotopes in surface snow and the overlying atmosphere. The deposited nitrate flux FD represents the state of nitrate that has just deposited onto the surface snow via dry deposition of gaseous nitrate or wet scavenge from the atmosphere and is close to the definition of the skin layer of snowpack, i.e., the uppermost several millimeters of surface snow (Erbland et al., 2013; Winton et al., 2020). Thus, if there are sufficient high-resolution skin layer observations, a direct comparison with the model output can be performed (i.e., FD vs. skin layer measurements). Moreover, since FD originates from the local atmosphere, if the air–snow nitrate transfer function (i.e., the mass and isotope relationships between atmospheric nitrate and the deposited nitrate) is known, the calculated FD could be used to infer the state of local atmospheric nitrate. In this study, the isotope transfer function is applied instead of the mass transfer function because of its simplicity, especially for Δ¹⁷O, which is assumed to be conserved during deposition owing to its mass-independent nature. For δ¹⁵N, we assume that the deposition of atmospheric nitrate is associated with a fractionation constant (ε_d) of +10 ‰ following Erbland et al. (2013). We can either directly compare the modeled isotopes of FD with the observed values in the skin layer or with local atmospheric signals by including the differences (only for δ¹⁵N) between FD and atmospheric nitrate.

In this study, we chose two typical polar sites, Summit, Greenland, and Dome C, Antarctica, to conduct case studies in order to test the performance of the inverse model. These two sites were chosen for several reasons. First, these two sites represent typical polar sites with both relatively high (Summit) and extremely low (Dome C) snow accumulation rates. Second, there are sufficient atmospheric and/or snow observations at these two sites, which informs model input parameters and allows for comparison of the model results with observations. Third, these two sites are hot spots of ice-core drilling, and future work using the inverse model on ice-core nitrate records from these sites can be performed. In addition, there have already been studies simulating the post-depositional processing of snow nitrate at these two sites by using the forward TRANSITS model (Erbland et al., 2015; Jiang et al., 2021). Most of the model parameters in this study are kept the same as the original TRANSITS simulation unless otherwise mentioned. The major parameters used in this study are summarized in Table 2. Below, we specifically describe how we chose the initial model values/conditions for simulations at these two sites.

Table 2 Values of major parameters used in the model simulations at two different sites.

^a Annual average value: the weekly resolution data were adopted from Jiang et al. (2022). ^b Adjusted according to the best fit of snowpack nitrate δ¹⁵N profile at Dome C (Appendix C).

3.1 Summit, Greenland

Summit, Greenland, is a typical site with high snow accumulation rate (250 kg m⁻² a⁻¹; Dibb and Fahnestock, 2004) at present, and weekly resolved snow accumulation data exists (Burkhart et al., 2004), allowing for the precise dating of the snowpack nitrate profile (Jiang et al., 2022). The snowpack nitrate concentration and isotope data with weekly resolution at Summit compiled in Jiang et al. (2022) were used as initial model input values to represent the archived snow nitrate signals.

3.2 Dome C, Antarctica

The present snow accumulation rate at Dome C, Antarctica, is extremely low (28 kg m⁻² a⁻¹, Erbland et al., 2013), and it is currently impossible to discern seasonal or sub-seasonal nitrate patterns owing to the limited resolution of snowpack measurements. Erbland et al. (2013) reported five snowpack nitrate depth profiles at Dome C that extended just below the photic zone. To predict the final archived nitrate concentration and isotopes, Erbland et al. (2013) fitted these depth profiles with an exponential function, and the obtained asymptotic values were regarded as the final preserved nitrate signal. The average asymptotic values for the five snowpacks were (21.2 ± 18.1) ng g⁻¹, (273.6 ± 64.0) ‰, and (26.0 ± 1.9) ‰ for nitrate concentration, δ¹⁵N and Δ¹⁷O, respectively. These values were used as the annual averages of the preserved nitrate at Dome C in this study.

We note the seasonality of the archived nitrate concentration is important because it determines the magnitude of FP and FD at each time step in the model. In simulations of Dome C, we designed three cases with different weekly concentration distributions in a year. In case 1, the weekly nitrate concentrations were assumed to be uniform throughout a year. In case 2, the weekly archival nitrate concentrations were assumed to be a Gaussian-type distribution to match the observed seasonality in skin layer nitrate concentrations at Dome C (Erbland et al., 2013):

$$\begin{array}{}\text{(20)}& c\left(n\right)=c_{\mathrm{a}}\times \left(a+b\times \exp \left(-{\displaystyle \frac{{\left(n-n_{\mathrm{0}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right)\right).\end{array}$$

In Eq. (20), c_a represents the annual average snow nitrate concentration, n represents the week number (1 to 52, here week 1 is defined as the first week in January for the Northern Hemisphere sites or the first week in July for the Southern Hemisphere sites), and the shape parameters (a, b, σ) were determined by the best fit of skin layer nitrate concentrations (Appendix D). n₀ represents the week when nitrate concentration peaks in a year and was set to be 26 according to the observed maximum nitrate concentrations in the skin layer in local midsummer (Erbland et al., 2013). However, since nitrate deposited in different weeks of a year would have experienced different amounts of total actinic flux and nitrate deposited in autumn undergoes minimal degree of photolysis (Jiang et al., 2022), it is likely that the summer peak would shift toward autumn by final preservation. As such, we also prescribed a “shifted peak” distribution in case 3, and in this case n₀ was set equal to 35 in Eq. (20), while other parameters were the same as in case 2.

To determine the uncertainties in the model results caused by these artificially assumed nitrate profiles, we applied a Monte Carlo method; i.e., the exact initial value in snow at each week was set arbitrarily as follows:

$$\begin{array}{}\text{(21)}& c_{\mathrm{r}}=c_{\mathrm{a}}+U\left(-\mathit{\sigma },\mathit{\sigma }\right),\end{array}$$

where c_a represents the prescribed initial value of annual-mean snow nitrate concentration in each case as described above, U represents a uniformly distributed random variable, and σ represents the standard error of the observed c_a. The obtained time series with random error was normalized again as final model inputs. All three cases were repeated 1000 times, and the model results were used to evaluate the uncertainties.

For isotopic ratios (δ¹⁵N and Δ¹⁷O) of the archived nitrate, their seasonality was omitted in this study to simplify the model calculations, and for the results at Dome C we only compared the modeled results with observations at annual scale given the unknown seasonal inputs of these parameters. Note that reconstruction of atmospheric signals of δ¹⁵N and Δ¹⁷O from ice-core records usually use a coarser resolution than sub-annual variations, justifying our annual averaging approach.

4 Results and discussion

4.1 Model results at Summit, Greenland

4.1.1 Comparison of local atmospheric variations

Currently there are no skin layer observations at Summit, so we used the monthly atmospheric nitrate isotopic data from aerosol observations at Summit reported by Jiang et al. (2022) to compare with the modeled atmospheric nitrate isotopic variations. To reduce the uncertainty of model results owing to uncertainties associated with the weekly dating of snowpack, we only compared the modeled monthly averages with observations, and the uncertainties of the monthly model results were calculated as 1 standard error of the mean of results from different weeks. As shown in Fig. 3, the modeled seasonality in atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) generally agrees well with the observed seasonal variations, while for δ¹⁵N(NO${}_{\mathrm{3}}^{-}$), the model predicted a similar seasonality as the observations, though in the winter half year the model underestimated the absolute values in comparison with the observations. In addition, the modeled and observed atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) are both close to snowpack Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$). This is as expected at Summit since the deposition of atmospheric nitrate is assumed to be a mass-dependent process. The only process that can alter Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) in snow is the cage effect (McCabe et al., 2005), which is negligible under present Summit conditions. The inverse model calculated a small cage effect of 0.15 ‰ on Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) by comparing the annual weighted average of FA and FD, which is close to the value of 0.19 ‰ predicted by the TRANSITS model (Jiang et al., 2021). At an annual scale, the modeled and observed average atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) are $-\mathrm{17.5}\pm \mathrm{3.0}$ ‰ and −14.8 ± 7.3 ‰, while for Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) the values are 28.8 ± 2.6 ‰ and 28.6 ± 3.2 ‰ respectively, suggesting that the inverse model reproduced the atmospheric observations quite well.

Figure 3 Comparison between the modeled (red dots) and observed (blue stars) seasonal variations in atmospheric nitrate (a) δ¹⁵N and (b) Δ¹⁷O at Summit, Greenland. The dashed line with gray triangle represents the snowpack measurements as model inputs (i.e., the monthly archived snow nitrate δ¹⁵N and Δ¹⁷O from snowpack observations) (Jiang et al., 2022). The atmospheric observations were from Jiang et al. (2022).

For δ¹⁵N(NO${}_{\mathrm{3}}^{-}$), the modeled atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) seasonality is comparable to the observations, but the absolute values display some discrepancies in autumn and winter. In particular, the modeled and observed average δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) values in the summer half year (from March to August) are −17.6 ± 3.5 ‰ and −16.0 ± 7.8 ‰, respectively, while in the winter half-year they are −16.0 ± 7.8 ‰ and −12.0 ± 4.1 ‰, respectively. The model–observation difference in the winter half year may be related to the model setup of a constant ε_d of +10 ‰. As discussed by Jiang et al. (2022), the partition between nitrate deposition mechanisms (i.e., wet vs. dry deposition) may result in seasonally different air–snow transfer functions for δ¹⁵N(NO${}_{\mathrm{3}}^{-}$). It has been observed that δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) of dry deposition is generally higher than wet deposition (Beyn et al., 2014; Heaton, 1987). This implies that dry deposition likely possesses a larger ε_d, perhaps because wet deposition can scavenge all or most of atmospheric nitrate, leading to small to no isotope fractionation. Given the potential seasonal changes in the relative fraction of dry versus wet deposition at Summit, using a constant ε_d in the model would likely cause discrepancies in one season but not in the other. Some observations at Summit indicate that snowfall activities are more frequent and severe in summer months (June–September) than in winter (Castellani et al., 2015; Bennartz et al., 2019), which implies that dry deposition of atmospheric nitrate is more important in winter instead of summer. Thus, the model–observation discrepancies in the winter half year cannot be explained by the seasonal shift in the ratio between wet and dry deposition, as more dry deposition may result in a larger isotope effect in winter.

Alternatively, we note the ε_d itself may have a seasonality which could be caused by the temperature dependence of nitrate absorption onto ice grains (Abbatt, 1997) or be influenced by other mechanisms such as the stability of the boundary layer. In fact, observations at Dome C indicated the averaged enrichment in skin layer δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) related to atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) is +25 ‰ in summer, while in winter the value is +10 ‰ (Erbland et al., 2013). This may indicate a larger ε_d in summer than in winter, though the summer skin layer δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) is probably more or less influenced by photolysis which tends to increase δ¹⁵N(NO${}_{\mathrm{3}}^{-}$). However, in the inverse model, ε_d was set as +10 ‰ throughout the year for Summit (note this value is consistent with the observed difference between surface snow and local atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) at Summit in May and June; Fibiger et al., 2016). If at Summit the ε_d in winter is lower than that in summer, the modeled average δ¹⁵N(NO${}_{\mathrm{3}}^{-})$ in the winter half year would have been underestimated. This at least explains in part the model–observation discrepancies in winter half year δ¹⁵N(NO${}_{\mathrm{3}}^{-}$). Future work on the degree of nitrogen isotope fractionation during atmospheric nitrate deposition and the causal factors is necessary to further investigate this issue.

4.1.2 Flux and isotopes of primary nitrate

The major function of the inverse model is to reconstruct the primary nitrate flux and its isotopes by using ice-core nitrate records. Primary nitrate flux is closely associated with atmospheric nitrate on the regional scale and could be further linked to the atmospheric abundance of its precursor NO_x. The isotopic composition of F_pri could provide extra information. For example, the δ¹⁵N of F_pri may be used to infer the variations in NO_x source emissions if other factors influencing isotope fractionation during the atmospheric conversion of NO_x to nitrate can be constrained. The Δ¹⁷O of F_pri depends on the relative concentration of major atmospheric oxidations such as O₃ and HO₂/RO₂ radicals and thus could be used to reflect regional atmospheric oxidation environment (Geng et al., 2017; Sofen et al., 2014).

Figure 4 The modeled monthly (a) δ¹⁵N, (b) Δ¹⁷O, and (c) fluxes of primary nitrate (F_pri) to Summit, Greenland. The concentration, δ¹⁵N, and Δ¹⁷O values of the archived snow nitrate as model inputs are also shown for comparison (gray triangle with dashed line).

The model-calculated seasonal variations in F_pri to Summit are shown in Fig. 4. The annual flux of primary nitrate was calculated to be 6.96 × 10⁻⁶ kgN m² a⁻¹, which is in the same range (≈ 2–3 × 10⁻⁶ kgN m² a⁻¹) as model results from Zatko et al. (2016) using the GEOS-Chem model and is about 2 orders of magnitude lower than the atmospheric nitrate deposition flux in the mid-latitude area (Gao et al., 2019; Zhang et al., 2012). The seasonality of F_pri displays a bimodal mode with a major summer peak and a secondary peak in late winter/early spring, in contrast the preserved snowpack nitrate concentration which peaks in spring/summer. The maximum F_pri in summer could be caused by the enhanced temperature-dependent precursor NO_x emissions such as from soil microbes (Pilegaard et al., 2006) as well as the more active photochemistry in summer, both of which would promote more efficient atmospheric nitrate production. It is interesting that the secondary F_pri peak in early spring is coincident with the timing of the spring Arctic haze phenomenon (Quinn et al., 2007), as well as the occasional spring nitrate concentration peak in snowpack and ice cores at Summit (Geng et al., 2014), though the exact timing of the seasonal peaks needs further investigation.

The modeled Δ¹⁷O of F_pri is close to the measurements in the snowpack with minimum values in summer, suggesting the Δ¹⁷O signal of primary nitrate is well preserved under current Summit conditions. The seasonal variations in Δ¹⁷O of F_pri can be understood in terms of the different production mechanisms of atmospheric nitrate (Alexander et al., 2020). In summer, ample solar radiation enhances the photochemical production of HNO₃ from the NO₂ + OH pathway, the Δ¹⁷O of which is lowest compared with other nitrate formation pathways. While in winter, the dominant N₂O₅ hydrolysis pathway produces nitrate with high Δ¹⁷O. Such seasonal patterns have been widely observed globally as summarized in Alexander et al. (2020).

The modeled δ¹⁵N of F_pri ranges from −10.3 ‰ to 5.0 ‰, which falls well within the reported atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) values in continental and marine boundary layer in both hemispheres in regions not impacted by snowpack emission (Li et al., 2022; Lim et al., 2022; Morin et al., 2009; Shi et al., 2021). However, the seasonal pattern of δ¹⁵N of F_pri, which displays a summer maximum, is opposite to the typical seasonal pattern of atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) found in mid-latitude continental areas, where higher δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) values in winter and lower δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) values in summer are widely observed (e.g., Beyn et al., 2014; Freyer, 1991; Fang et al., 2021; Lim et al., 2022; Esquivel Hernández et al., 2023). This summer high and winter low δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) in F_pri is instead consistent with the observations at two Arctic coastal sites (Morin et al., 2012, 2008), where the summer high atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) is strongly correlated with air temperature. Morin et al. (2008) suggest the δ¹⁵N(NO${}_{\mathrm{3}}^{-}$)–temperature relationship observed at the Arctic coastal sites may be related to physicochemical transformations of nitrate in the Arctic and during the transport of nitrate and reactive nitrogen from the mid-latitudes, though the specific mechanism is unknown.

Another possibility to explain the higher modeled summer δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) in F_pri is that there may be more anthropogenic nitrate transported from mid-latitudes to Greenland in summer than in winter. F_pri is comprised of nitrate originating from the mid-latitudes as well as nitrate formed in the Arctic region. Morin et al. (2009) suggested that air parcels originating from regions with more anthropogenic impacts carries nitrate with higher δ¹⁵N, which was confirmed by subsequent studies (Li et al., 2022; Vicars and Savarino, 2014; Shi et al., 2021). The increased frequency of air sources originating from the North America in summer compared to winter (Kahl et al., 1997) could thus lead to more anthropogenic nitrate to Greenland in summer, resulting in higher summer δ¹⁵N of primary nitrate than winter.

The potential link between δ¹⁵N of F_pri and its precursor NO_x emissions is not further discussed here, as recent studies have shown that the isotopic effect during NO_x photo-recycling is complex (Li et al., 2020) and may dominate δ¹⁵N variations in atmospheric nitrate (Fang et al., 2021; Li et al., 2021). More comprehensive studies on the isotopic effects during atmospheric nitrate formation as well as the potential fractionation during transport are required to further link δ¹⁵N of F_pri with its precursors and/or source regions.

4.2 Model results at Dome C, Antarctica

4.2.1 Snowpack nitrate profile in the photic zone

Since Dome C snowpack exhibits very distinct trends in the concentration and isotopic ratio of nitrate in the photic zone, we first examine the modeled summer snowpack nitrate profile at Dome C in comparison with the previous observations (Erbland et al., 2013; Frey et al., 2009) in Fig. 5. The TRANSITS model results (Erbland et al., 2015) are also shown in Fig. 5 for comparison. Both models reproduce the observed decrease in nitrate concentration and the large enrichments in δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) well. We note that the predicted surface snow nitrate concentration is higher than the observations by both models. This is because the modeled concentration represents the state that atmospheric nitrate has just deposited onto the snow surface, while the observed skin layer snow may have already undergone snow metamorphism and/or post-depositional processing (Winton et al., 2020). This is also supported by recent observations at Dome C that newly deposited diamond dust could possess nitrate concentrations up to 2000 ppb (Winton et al., 2020), within the range of model predictions.

Figure 5 Comparison between the observed and modeled snowpack nitrate concentrations, δ¹⁵N, and Δ¹⁷O at Dome C. The red lines with circles represent four observed snowpack nitrate profiles at Dome C from Erbland et al. (2013) and Frey et al. (2009), while the blue and black lines are modeled results from the forward model (i.e., the TRANSITS model) and the inverse model in this study, respectively. The yellow background represents the depth of the photic zone.

The decreasing trend in Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) within the photic zone is also reproduced by these two models, caused by the cage effect during nitrate photolysis. However, the TRANSITS model appears to underestimate snowpack Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$), while the inverse model performs better in snowpack Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) simulation. This is because in the TRANSITS model, snow Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) is controlled by a combination of Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FD and the subsequent cage effect after deposition. At Dome C, Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FD is dominated by locally formed atmospheric nitrate (i.e., FP) (Erbland et al., 2015), which is in turn determined by the prescribed Δ¹⁷O transfer during NO–NO₂ cycling and the subsequent OH oxidation of NO₂ under sunlight conditions in the model. However, Savarino et al. (2016) demonstrated that the standard chemistry (i.e., exclusive oxidation of NO₂ by OH in summer) and the associated isotopic mass balance applied to Δ¹⁷O (i.e., the one used by the direct model) does not hold in the Dome C atmosphere, with this standard approach systematically underestimating the observations. Our inverse model is in line with this conclusion. The inverse model calculates atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) from the archived snow Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) by subtracting the cage effect but does not assume any specific chemical reaction in the atmospheric box, contrary to the forward model. Therefore, although the inverse model uses the same method as the TRANSITS model to calculate Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of locally formed atmospheric nitrate (FP), it does not include any hypothesis of how local nitrate is formed. The good match between observations and inverse model output is a further demonstration that atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) is not in agreement with the standard daylight chemistry of nitrate formation (i.e., NO₂ + OH).

Overall, the consistency of the modeled and observed snowpack nitrate profiles suggests that the effect of post-depositional processing is properly represented by the inverse model. This confirms that the inverse model can properly reproduce snow nitrate concentrations and isotopes in the photic zone, which are intermediate statuses between archived and atmospheric nitrate.

4.2.2 Skin layer and atmospheric nitrate

Although the time step in the model was set to 1 week for Dome C, we mainly focus on the annual average values since the seasonal information of the archived nitrate profiles is unknown. We note that the modeled isotopic compositions of snowpack and skin layer nitrate are not affected by the prescribed nitrate concentration seasonality. This is because the total nitrate loss fraction and the induced isotopic effect only depend on the total amount of actinic flux received during snow burial. In the following discussion we only report and discuss the modeled isotopes of local atmospheric and skin layer nitrate from case 1; i.e., the archived snow nitrate concentration was assumed to be constant throughout the year.

Figure 6 Comparison between the observed and modeled annual averages of δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) in the atmosphere and snow skin layer at Dome C. The solid line in the box plot indicates the median value, while the dash line represents the average value.

The observed annual average δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) values in the skin layer at Dome C are 18.0 ± 11.7 ‰ and 33.6 ± 1.4 ‰, respectively (Erbland et al., 2013), while the modeled skin-layer values are 15.7 ± 38.6 ‰ and 33.3 ± 4.7 ‰, respectively, in good agreement with the observations. The observed annual average atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) values are −1.3 ± 11.6 ‰ and 31.4 ± 4.6 ‰, while the modeled values are 8.0 ± 11.7 ‰ and 33.6 ± 1.4 ‰, respectively. Note the average observed δ¹⁵N values in this study were calculated as arithmetic mean instead of mass-weighted mean reported in Erbland et al. (2013) since the inverse model cannot directly calculate the nitrate concentration in the atmosphere. Nevertheless, the modeled averages are similar to the observed averages except for atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$). The difference between the modeled and observed atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) could be again related to the constant ε_d used in the model. As discussed earlier, in the model we followed Erbland et al. (2015) to set ε_d=10 ‰ throughout the year, while observations at Dome C indicate that ε_d could be as large as 25 ‰ in summer instead of 10 ‰ (Erbland et al., 2013). Hence the modeled atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) could be overestimated. This reinforces that it is necessary to further explore the isotope effects on δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) during atmospheric nitrate deposition.

4.2.3 Flux and isotopes of primary nitrate

In Fig. 7, the flux of primary nitrate (F_pri) and its mean isotopes from the three difference cases (i.e., different nitrate concentration seasonality in archived snow) are displayed. Similar to the previous section, we only focus on their annual means. Note when calculating Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri, Δ¹⁷O values of locally formed nitrate are necessary. As discussed in Sect. 4.2.1, following the method in the TRANSITS model to calculate Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of locally formed atmospheric nitrate would underestimate Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FP at Dome C. This is especially evident in summer when the snow sourced nitrate (i.e., FP) dominates the atmospheric nitrate budget, and the calculated Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FP by the TRANSITS model is about 6 ‰ lower than observed atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) (Erbland et al., 2015). Thus, in the inverse model, when calculating Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri at Dome C, Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FP was not calculated using the model default method as in the TRANSITS but prescribed as the observed atmospheric summer Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$). Otherwise, the modeled Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of F_pri would be higher than 70 ‰, which is highly unrealistic. Note that this is not an issue at Summit, Greenland, because FP does not dominate the atmospheric nitrate budget in summer there.

Figure 7 Frequency histogram of the calculated primary nitrate flux and its mean δ¹⁵N/Δ¹⁷O values at Dome C under three different archival snow nitrate concentration distributions: (a–c) Case 1 – uniform distribution, (d–f) Case 2 – Gaussian-type function, and (g–i) Case 3 – shifted Gaussian-type function. The solid black lines represent the fitted Gaussian function of the frequency distribution. The dashed blue lines represent the mean values of F_pri and its δ¹⁵N and Δ¹⁷O, which are labeled on the top of each subplot.

As shown in Fig. 7, although the prescribed archived nitrate concentration seasonality does not alter the modeled snowpack and atmospheric nitrate isotopes, it has a profound impact on the modeled primary nitrate flux and its isotopes. In particular, under the three cases of different seasonal distributions of the archived snow nitrate concentrations, the modeled F_pri and its annual mean δ¹⁵N and Δ¹⁷O range from 1.5 to 2.2 × 10⁻⁶ kgN m⁻² a⁻¹, 6.2 ‰ to 29.3 ‰, and 48.8 ‰ to 52.6 ‰, respectively. The inverse model-calculated F_pri is smaller than the value used in the original TRANSITS (8.2 × 10⁻⁶ kgN m⁻² a⁻¹) in Erbland et al. (2015), but this is easily resolved given the large uncertainty in the archived nitrate concentration used as model input.

The modeled annual mean δ¹⁵N of F_pri ranges from 6.2 ‰–29.3 ‰, in contrast with the observed atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) in southern mid-latitude area or the Southern Ocean where δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) is in general negative or close to 0 (Morin et al., 2009; Shi et al., 2018, 2021). The modeled positive δ¹⁵N of F_pri is however consistent with the wintertime atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) observed in Antarctica when the effect of photolysis is null and local atmospheric nitrate likely reflects F_pri. The maximum atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) in winter was found to be 10.8 ‰ at Dumont d’Urville station (Savarino et al., 2007), 12.8 ‰ at Dome C (Erbland et al., 2013), and 13.9 ‰ at Zhongshan station (Shi et al., 2022). These positive δ¹⁵N values have been linked to stratospheric denitrification as nitrate produced in stratosphere is suggested to be 19 ± 3 ‰ by considering the fractionation induced by different of N₂O photolysis channels (Savarino et al., 2007). Therefore, the modeled flux and δ¹⁵N of F_pri point towards the dominance of stratospheric denitrification in nitrate budget at Dome C.

The modeled Δ¹⁷O of F_pri is also very high for all three cases (48.8 ‰–52.6 ‰). The measured bulk Δ¹⁷O of surface ozone in Antarctica is about 26 ‰ (Ishino et al., 2017; Savarino et al., 2016) that fits well with the global tropospheric average of 25.4 ‰ (Vicars and Savarino, 2014). Given that the oxygen mass-independent fractionation signal of ozone is mainly occupied by the terminal oxygen atom and transferred to other molecular, atmospheric nitrate of tropospheric origin should possess a Δ¹⁷O signal less than 39 ‰ (Mauersberger et al., 2003; Savarino et al., 2008), which cannot explain our calculated high Δ¹⁷O of F_pri. However, the bulk Δ¹⁷O of stratospheric ozone was measured to be 34.3 ± 3.6 ‰ (Lämmerzahl et al., 2002; Krankowsky et al., 2000), which indicated that nitrate produced in the stratosphere could gain a higher Δ¹⁷O signature from ozone (Lyons, 2001). It has been observed in Antarctica that the atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) could exceed 40 ‰ in winter and early spring when stratospheric denitrification occurs (Ishino et al., 2017; Walters et al., 2019; Erbland et al., 2013; Savarino et al., 2007; Shi et al., 2022). A recent study also revealed that the surface snow Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) at Dome C frequently exceeds 40 ‰ during winter/spring and could sometimes reach up to 50 ‰ (Akers et al., 2022a). As we mentioned previously, these winter Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) observations likely reflect the primary nitrate signal at that time since the photolysis of snow nitrate does not occur due to lack of sunlight. Thus, the high modeled Δ¹⁷O of F_pri seems to again indicate a dominant role of stratosphere denitrification in external nitrate source to Dome C, similar to what can be reflected from the modeled δ¹⁵N of primary nitrate. In addition, Erbland et al. (2015) estimated that stratospheric denitrification nitrate flux is (4.1 ± 2.5) × 10⁻⁶ kgN m⁻² a⁻¹ in Antarctica, while our calculated F_pri of 1.5–2.2 × 10⁻⁶ kgN m⁻² a⁻¹ at Dome C is within the same range.

In sum, we acknowledge that there are many factors that would affect the model results, such as the archived snow nitrate concentration and isotopes, the export fraction (f_exp), and the cage effect fraction (f_c). These need to be further explored by observations to improve the model performance.

5 Model sensitivity tests: the impact of f_exp and f_c

In this section, we report the sensitivity test results to elucidate the impact of two model parameters that lack direct observational constraints: the export fraction (f_exp) and the cage effect fraction (f_c). We mainly focus on the annual net loss and the differences in isotopes between F_pri and FA in accordance with the resolution of ice-core measurements.

Figure 8 Model sensitivity test results of two parameters f_exp and f_c for Summit (a–c) and Dome C (d–f). The annual nitrate loss fraction is defined as $\mathrm{1}-\mathrm{FA}/F_{\mathrm{pri}}$ following Jiang et al. (2021). Red stars represent the values of f_exp and f_c used in model simulations.

The sensitivity test results are shown in Fig. 8. The annual loss fraction (defined as $\mathrm{1}-\mathrm{FA}/F_{\mathrm{pri}}$) represents the final preservation of primary nitrate after post-depositional processing (Jiang et al., 2021). The inverse model predicts an annual loss fraction of 3.5 % under present Summit conditions, which is close to the TRANSITS model prediction of 4.1 % (Jiang et al., 2021). This small discrepancy is likely caused by the use of simplified snow radiative transfer parameterization in the inverse model. In addition, the differences of δ¹⁵N and Δ¹⁷O between FA and F_pri are also in good agreement with the TRANSITS model. As expected, larger f_exp and smaller f_c would result in a higher degree of net loss in F_pri and larger isotopic effects in both δ¹⁵N and Δ¹⁷O. However, under present-day conditions at Summit, the preserved snow nitrate concentrations and isotopes at the annual scale are only altered slightly and the degree of changes is insensitive to f_exp and f_c.

For Dome C, the model results are sensitive to f_exp when f_c is small and become sensitive to f_c when f_exp is larger. In particular, the isotopic signature is largely independent of f_exp when f_exp is larger than 0.1–0.2. In addition, the Δ¹⁷O results display a non-monotonic response to these two parameters, especially when f_exp approaches zero (Fig. 8f). A similar phenomenon was seen in the TRANSITS model simulations in Erbland et al. (2015), where they found that the model results could not converge when f_exp was set to zero. The high sensitivity of model parameters renders difficult to reconstruct the historical variations in primary nitrate based on ice-core records at Dome C unless these parameters are precisely constrained. For present-day conditions, f_exp and f_c could be constrained by atmospheric and snowpack observations (Erbland et al., 2015), but it is unknown whether these values could be applied to different climate conditions. In addition, the difficulties in choosing an appropriate archival nitrate concentration profile as model initial conditions would add extra uncertainties to the model results.

6 Conclusions and implications

In this study, we introduce an inverse model which is designed to correct for the effects of post-depositional processing on ice-core nitrate concentration and its isotopes. The model was tested against present-day Summit, Greenland, and Dome C, Antarctica, conditions to validate its performance under different snow accumulation rates. Model results compared to observations demonstrate that the inverse model is capable of adequately correcting the effect of post-depositional processing. The modeled atmospheric nitrate δ¹⁵N/Δ¹⁷O values at Summit are generally in good agreement with observations but with a slight underestimate in winter δ¹⁵N(NO${}_{\mathrm{3}}^{-}$), which is likely because the model does not treat correctly the seasonal differences in nitrogen isotope fractionation during deposition (ε_d). At Dome C, the model also well reproduced the observed snowpack nitrate profiles in the photic zone, the annual skin layer δ¹⁵N/Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$), and atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) at Dome C but again overestimated the average atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) probably also due a low bias in ε_d used in the model. A better quantification on the isotope fractionation of δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) during deposition is therefore needed.

The inverse working flow of this new model also enables us to qualitatively retrieve information regarding primary nitrate from the archived snow nitrate. The calculated seasonality in δ¹⁵N of F_pri at Summit displays a maximum in mid-summer that is distinct from the observed spring δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) peak in snowpack. This seasonal pattern is in contrast with observed atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) variations in mid-latitudes, which is thought to be the major aerosol source region to Summit, but is consistent with the atmospheric δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) variations observed in the high-latitude Arctic region. The δ¹⁵N of F_pri may reflect seasonally varied main source regions to Greenland or a dominate role of high-latitude nitrate transport to Summit. At Dome C, both the magnitude of F_pri and its δ¹⁵N/Δ¹⁷O indicate a dominant role of stratospheric denitrification on nitrate budget at Dome C.

The inverse model is designed to help interpret ice-core nitrate records. Applying the inverse model to ice-core nitrate records needs knowledge of initial conditions. In particular, archived snow nitrate concentration and its δ¹⁵N and Δ¹⁷O, the snow accumulation rate, and light absorption impurity concentrations should be known for a given ice core. In addition, chemistry–climate models such as the ICECAP or GCAP model (Murray et al., 2021, 2014) would be also necessary to provide extra constraints, such as the oxidizing agent concentrations, total column ozone (TCO), wind field, and boundary layer heights for the past climates and are required to estimate Δ¹⁷O of FP and f_exp (Alexander et al., 2020; Jiang et al., 2021). The calculated primary nitrate flux and its δ¹⁵N and Δ¹⁷O can be further combined with the chemistry–climate model results to interpret its climate implications such as the variations in tropospheric NO_x and oxidant abundance, which would improve our understanding of key factors controlling the variability in atmospheric oxidation capacity under different climates.

Appendix A: Derivation of the nitrate mass and isotopic balance equations during photolysis (Eqs. 7–9)

In the inverse model, we follow Erbland et al. (2015) to separate the photolysis process of nitrate on an ice grain into two steps, i.e., the direct photolysis followed by the subsequent cage effect (Fig. A1). It is well-documented that secondary chemistry can occur during snow nitrate photolysis to reform nitrate and alter the isotopes of the remaining nitrate, which is termed as the cage effect (McCabe et al., 2005; Meusinger et al., 2014). To quantify this effect, Erbland et al. (2015) assigned an empirical parameter (f_c) to represent the fraction of the nitrate photoproducts that would undergo the cage effect and derived a value of 0.15 for f_c based on the observed decreasing trend of Δ¹⁷O (NO${}_{\mathrm{3}}^{-}$) in the snowpack at Dome C. As shown in Fig. A1, assuming a fraction (f_p) of initial snow nitrate was photolyzed and a fraction (f_c) of these photolyzed nitrate undergone cage effect, the mass balance equation for snow nitrate can be written as

$$\begin{array}{}\text{(A1)}& c\left(\mathrm{SN}\right)=c\left({\mathrm{SN}}^{\prime }\right)\left(\mathrm{1}-f_{\mathrm{p}}+f_{\mathrm{c}}{f}_{\mathrm{p}}\right),\end{array}$$

where the superscript represents the state before photolysis. It can be easily seen that Eq. (A1) is equal to Eq. (7).

The impact of the direct photolysis of snow nitrate on δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) can be described by the Rayleigh equation. We define the first-order photolysis rate constant of ¹⁴NO${}_{\mathrm{3}}^{-}$ and ¹⁵NO${}_{\mathrm{3}}^{-}$ as J and J^∗ and their concentration in snow as c and c^∗ respectively. The chemical kinetic equations of c and c^∗ can be represented as follows:

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\mathrm{d}c}{\mathrm{d}t}}=-Jc\text{(A3)}& {\displaystyle \frac{\mathrm{d}{c}^{\ast }}{\mathrm{d}t}}=-J^{\ast }{c}^{\ast }.\end{array}$$

Integrating Eqs. (A2) and (A3) yields Eqs. (A4) and (A5):

$$\begin{array}{}\text{(A4)}& {\displaystyle }c\left(t\right)=c\left(\mathrm{0}\right)e^{-{\int }_{\mathrm{0}}^{t}J\mathrm{d}t}\text{(A5)}& {\displaystyle }{c}^{\ast }\left(t\right)=c^{\ast }\left(\mathrm{0}\right)e^{-{\int }_{\mathrm{0}}^t{J}^{\ast }\mathrm{d}t}.\end{array}$$

Here c(0) represents the initial concentration before photolysis. The evolution of the isotopic ratio R, which is defined as the ratio of c and c^∗, follows Eq. (A6):

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}R\left(t\right)& ={\displaystyle \frac{c^{\ast }\left(t\right)}{c\left(t\right)}}={\displaystyle \frac{c^{\ast }\left(\mathrm{0}\right)}{c\left(\mathrm{0}\right)}}{e}^{-{\int }_{\mathrm{0}}^t(J^{\ast }-J)\mathrm{d}t}\\ & =R\left(\mathrm{0}\right)e^{-{\int }_{\mathrm{0}}^t(J^{\ast }-J)\mathrm{d}t}.\end{array}\end{array}$$

Since the delta value δ¹⁵N equals $R_{\mathrm{spl}}/R_{\mathrm{ref}}-\mathrm{1}$, where R_spl and R_ref refer to the isotope ratio of sample and standard respectively, Eq. (A6) can be further expanded to

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}\ln {\displaystyle \frac{\mathrm{1}+\mathit{\delta }\left(t\right)}{\mathrm{1}+\mathit{\delta }\left(\mathrm{0}\right)}}& =\ln {\displaystyle \frac{R\left(t\right)}{R\left(\mathrm{0}\right)}}=-\underset{\mathrm{0}}{\overset{t}{\int }}(J^{\ast }-J)\mathrm{d}t\\ & =-\underset{\mathrm{0}}{\overset{t}{\int }}J{\mathit{\epsilon }}_{\mathrm{p}}\mathrm{d}t=-\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}\underset{\mathrm{0}}{\overset{t}{\int }}J\mathrm{d}t\\ & =\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}\ln (\mathrm{1}-f_{\mathrm{p}}),\end{array}\end{array}$$

which is consistent with the form of the Rayleigh equation.

By applying the first-order Taylor expansion of $\ln (\mathrm{1}+{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}({\mathrm{NO}}_{\mathrm{3}}^{-}\left)\right)\approx {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}({\mathrm{NO}}_{\mathrm{3}}^{-}$), we obtain the relationship between the δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) before and after photolysis:

$$\begin{array}{}\text{(A8)}& \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_r\right)\approx \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}^{\prime }\right)-\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}\ln (\mathrm{1}-f_{\mathrm{p}}).\end{array}$$

The δ¹⁵N of the emitted NO₂ can be calculated via the mass balance equation:

$$\begin{array}{}\text{(A9)}& \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}^{\prime }\right)=\left(\mathrm{1}-f_{\mathrm{p}}\right)\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_r\right)+f_{\mathrm{p}}\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{NO}}_{\mathrm{2}}\right).\end{array}$$

Combining Eqs. (A8) and (A9) would yield

$$\begin{array}{}\text{(A10)}& \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{NO}}_{\mathrm{2}}\right)\approx \mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}^{\prime }\right)+{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}(\mathrm{1}-f_{\mathrm{p}})\ln (\mathrm{1}-f_{\mathrm{p}})}{f_{\mathrm{p}}}}.\end{array}$$

Because part of the photoproduct would undergo cage effect to reform nitrate (Fig. A1), the final state of snow δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) after photolysis can be calculated via isotopic mass balance equation:

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left(\mathrm{SN}\right)& ={\displaystyle \frac{\left(\mathrm{1}-f_{\mathrm{p}}\right)\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}_r\right)+f_{\mathrm{c}}{f}_{\mathrm{p}}\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{NO}}_{\mathrm{2}}\right)}{\mathrm{1}-f_{\mathrm{p}}+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\\ & =\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}\left({\mathrm{SN}}^{\prime }\right)-{\displaystyle \frac{\left(\mathrm{1}-f_{\mathrm{p}}\right)\left(\mathrm{1}-f_{\mathrm{c}}\right)\stackrel{\mathrm{‾}}{{\mathit{\epsilon }}_{\mathrm{p}}}\ln \left(\mathrm{1}-f_{\mathrm{p}}\right)}{\left(\mathrm{1}-f_{\mathrm{p}}\right)+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\end{array},\end{array}$$

which is equal to Eq. (8).

For Δ¹⁷O, it is assumed that direct photolysis of nitrate will only induce mass-dependent fractionation and has no impact on Δ¹⁷O. However, the produced NO₂ which is re-oxidized by an OH radical would lower its Δ¹⁷O by $\mathrm{2}/\mathrm{3}$ since an OH radical would rapidly achieve isotopic equilibrium with the surrounding water molecule and erases its Δ¹⁷O signal. Hence, by again using the isotopic mass balance equation, it can be shown that

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left(\mathrm{SN}\right)& ={\displaystyle \frac{\left(\mathrm{1}-f_{\mathrm{p}}\right)\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}^{\prime }\right)+\frac{\mathrm{2}}{\mathrm{3}}{f}_{\mathrm{c}}{f}_{\mathrm{p}}\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}^{\prime }\right)}{\mathrm{1}-f_{\mathrm{p}}+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\\ & ={\displaystyle \frac{\mathrm{1}-f_{\mathrm{p}}+\frac{\mathrm{2}}{\mathrm{3}}{f}_{\mathrm{c}}{f}_{\mathrm{p}}}{\mathrm{1}-f_{\mathrm{p}}+f_{\mathrm{c}}{f}_{\mathrm{p}}}}\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{SN}}^{\prime }\right),\end{array}\end{array}$$

which is equal to Eq. (9).

Figure A1 Sketch of the mass and isotopic transfer relationship during nitrate photolysis in snow grains. The italic black variables near the arrows indicate the fractional change in each subprocess.

Appendix B: Method for calculating Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FP

The Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) of FP is required to solve the mass balance equations. We follow the algorithm for calculating atmospheric Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) that has been widely used in a previous study (Alexander et al., 2020). Atmospheric NO₂ is assumed to rapidly achieve photochemical steady state (PSS) so that its Δ¹⁷O can be represented as follows:

$$\begin{array}{}\text{(B1)}& \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{NO}}_{\mathrm{2}}\right)={\displaystyle \frac{\begin{array}{c}{k}_{\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{O}}_{\mathrm{3}}\right]\\ +k_{\mathrm{NO}+\mathrm{BrO}}\left[\mathrm{BrO}\right]\end{array}}{\begin{array}{c}{k}_{\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{O}}_{\mathrm{3}}\right]+k_{\mathrm{NO}+\mathrm{BrO}}\left[\mathrm{BrO}\right]\\ +k_{\mathrm{NO}+{\mathrm{RO}}_{\mathrm{2}}}\left[{\mathrm{RO}}_{\mathrm{2}}\right]\end{array}}}\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{O}}_{\mathrm{3}}^{\ast }\right),\end{array}$$

where k represents different oxidation channels for NO in the atmosphere and Δ¹⁷O(O${}_{\mathrm{3}}^{\ast }$) represents Δ¹⁷O of the terminal oxygen in the ozone molecule, and RO₂ includes both HO₂ and other organic peroxyl radicals. Thus, to obtain Δ¹⁷O(NO₂), concentrations of ozone and oxidizing radicals are necessary. For the subsequent oxidation of NO₂, only the NO₂ + OH channel is considered, and Δ¹⁷O of OH is assumed to be zero owing to its isotopic equilibrium with atmospheric H₂O. Thus, Δ¹⁷O of the locally formed atmospheric nitrate (i.e., FP) can be calculated by

$$\begin{array}{}\text{(B2)}& \mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{NO}}_{\mathrm{3}}^{-}\right)={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\mathrm{\Delta }{}^{\mathrm{17}}\mathrm{O}\left({\mathrm{NO}}_{\mathrm{2}}\right).\end{array}$$

Appendix C: Adjusting the photolysis quantum yield used in Dome C simulations

The photolysis quantum yield (Φ) previously used in TRANSITS model simulation at Dome C was set to 0.026 in Erbland et al. (2015), which was adjusted to 0.015 in this study. Both values are obtained by best fitting the observed nitrate concentration and its δ¹⁵N profiles in the snowpack and are within the range of measurement results (0.003–0.44) for Dome C snow. However, using a value of 0.026 in the inverse model would severely overestimate the photolytic loss of snow nitrate, resulting in unrealistically high nitrate concentration in the skin layer (> 15 000 ppb) and excessive fractionation in δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) (500 ‰) and Δ¹⁷O(NO${}_{\mathrm{3}}^{-}$) (−18 ‰). Adjusting the quantum yield to 0.015 could well reproduce the observed nitrate concentration and δ¹⁵N(NO${}_{\mathrm{3}}^{-}$) in skin layer and snowpack by the inverse model.

The discrepancy in the chosen quantum yield in these two models is caused by whether diffusion of snow nitrate is included, as the diffusion process would tend to smooth the entire snowpack nitrate profiles and decrease the asymptotic values. The omission of the diffusion process in the inverse model is based on the following considerations. First, the snowpack nitrate profile at sites with even lower accumulation rates (Dome A in East Antarctica, with annual snow accumulation rate of about 19 kg m⁻² a⁻¹) does not display detectable smoothing on snowpack nitrate or its δ¹⁵N and Δ¹⁷O (Shi et al., 2015), suggesting that diffusion is not as important as assumed in the TRANSITS model. Second, the TRANSITS model would induce numerical diffusion during the division of snow layers in each time step, which results in rapid decreases in the amplitude of the simulated seasonality in snow, as can be seen in the simulated snow profiles in Winton et al. (2020).

Appendix D: Fitting the shape parameters for skin layer nitrate concentration at Dome C

To obtain a hypothetical seasonal pattern in the archived nitrate concentration, we assume that the normalized archived nitrate concentrations (by its arithmetic mean) would follow a same Gaussian-type distribution as nitrate concentrations in the skin layer, and the shape parameters (a, b, σ) are determined by using the least square regression method.

Figure D1 (a) Observed annual variations of skin layer nitrate concentration at Dome C (Erbland et al., 2013). (b) Normalized skin layer nitrate concentration by its arithmetic mean and the fitted Gaussian curve, with the shape parameters (a, b, σ) labeled in blue text.