Derivation of the reduced reaction mechanisms of ozone depletion events in the Arctic spring by using concentration sensitivity analysis and principal component analysis

The ozone depletion events (ODEs) in the springtime Arctic have been investigated since the 1980s. It is found that the depletion of ozone is highly associated with an auto-catalytic reaction cycle, which involves mostly the bromine-containing compounds. Moreover, bromide stored in various substrates in the Arctic such as the underlying surface covered by ice and snow can be also activated by the absorbed HOBr. Subsequently, this leads to an explosive increase of the bromine amount in the troposphere, which is called the “bromine explosion mechanism”. In the present study, a reaction scheme representing the chemistry of ozone depletion and halogen release is processed with two different mechanism reduction approaches, namely, the concentration sensitivity analysis and the principal component analysis. In the concentration sensitivity analysis, the interdependence of the mixing ratios of ozone and principal bromine species on the rate of each reaction in the ODE mechanism is identified. Furthermore, the most influential reactions in different time periods of ODEs are also revealed. By removing 11 reactions with the maximum absolute values of sensitivities lower than 10 %, a reduced reaction mechanism of ODEs is derived. The onsets of each time period of ODEs in simulations using the original reaction mechanism and the reduced reaction mechanism are identical while the maximum deviation of the mixing ratio of principal bromine species between different mechanisms is found to be less than 1 %. By performing the principal component analysis on an array of the sensitivity matrices, the dependence of a particular species concentration on a combination of the reaction rates in the mechanism is revealed. Redundant reactions are indicated by principal components corresponding to small eigenvalues and insignificant elements in principal components with large eigenvalues. Through this investigation, aside from the 11 reactions identified as unimportant in the concentration sensitivity analysis, additionally nine reactions were indicated to contribute only little to the total response of the system. Thus, they can be eliminated from the original reaction scheme. The results computed by applying the reduced reaction mechanism derived after the principal component analysis agree well with those by using the original reaction scheme. The maximum deviation of the mixing ratio of principal bromine species is found to be less than 10 %, which is guaranteed by the selection criterion adopted in the simplification process. Moreover, it is shown in the principal component analysis that O(1D) in the mechanism of ODEs is in quasi-steady state, which enables a following simplification of the reduced reaction mechanism obtained in the present study.

In the present study, a reaction scheme representing the chemistry of ozone depletion and halogen release is processed with two different mechanism reduction approaches, namely the concentration sensitivity analysis and the principal component analysis. In the concentration sensitivity analysis, the interdependence of the mixing ratios of ozone and principal bromine species on the rate of each reaction in the mechanism of ODEs is identified. Furthermore, the most influential reactions in different time periods of ODEs are also revealed. By removing 11 reactions with the maximum absolute values of sensitivities 10 lower than 10%, a reduced reaction mechanism of ODEs is derived. The onsets of each time period of ODEs in simulations using the original reaction mechanism and the reduced reaction mechanism are identical while the maximum deviation of the mixing ratio of principal bromine species between different mechanisms is found less than 1%.
By performing the principal component analysis on an array of the sensitivity matrices, the dependence of a particular species concentration on a combination of the reaction rates in the mechanism is revealed. Redundant reactions are indicated 15 by principal components corresponding to small eigenvalues and insignificant elements in principal components with large eigenvalues. Through this investigation, aside from the 11 reactions which have been identified as unimportant in the concentration sensitivity analysis, additionally nine reactions were identified to contribute only little to the total response of the system. Thus, they can be eliminated from the original reaction scheme. The results computed by applying the reduced reaction mechanism derived after the principal component analysis agree well with those by using the original reaction scheme. The 20 maximum deviation of the mixing ratio of principal bromine species is found less than 10%, which is guaranteed by the selection criterion adopted in the simplification process. Moreover, it is shown in the principal component analysis that O( 1 D) in

Introduction
Since the first discovery by Oltmans (1981) in observations conducted at Barrow, Alaska in 1977, the ozone depletion events (ODEs) in polar spring have attracted great attention. It was found by Oltmans (1981) that the surface ozone mixing ratio in 5 Barrow dropped from a background value (∼40 ppb, ppb = parts per billion) to a level lower than 1 ppb within a couple of days.
In a following measurement performed at Alert, Canada, not only the occurrence of the ozone depletion but also a negative correlation between the ozone mixing ratio and the concentration of filterable bromine (f-Br) was confirmed (Bottenheim et al., 1986;Barrie et al., 1988). Since then, ODEs and the associated bromine enhancement have been reported from various observation sites in polar regions (Kreher et al., 1997;Frieß et al., 2004;Jones et al., 2006;Wagner et al., 2007;Helmig et al., 10 2007; Jones et al., 2009Jones et al., , 2010Helmig et al., 2012;Halfacre et al., 2013).
Although, the concentration increase of the brominated species was found during ODEs, the mechanism responsible for the destruction of ozone with the involvement of bromine in the troposphere of Arctic remained a matter of debate until the role of bromine monoxide (BrO) was uncovered by Hausmann and Platt (1994). It was suggested by Hausmann and Platt (1994) that the bromine species participate in an auto-catalytic reaction cycle (I), namely through which the ozone in the boundary layer is consumed without any loss of bromine. In the reaction sequence (I), BrO participates in self reactions through which Br atoms are formed. Since the reaction between ozone and Br is rapid, ozone in the boundary layer is continuously consumed in the presence of sunlight via reaction sequence (I). However, the reaction cycle (I) is unable to explain the fast enhancement of bromine in the ambient air. Therefore, another reaction cycle (II) is proposed 20 as follows, In sequence (II), the HOBr, which is formed through the oxidation of BrO, is capable of activating bromide ions stored in various polar substrates such as the suspended aerosols and ice/snow-covered surfaces. As a result, the bromine amount in C 2 H 4 are critical for the ozone budget in the troposphere. Later on, a new comprehensive atmospheric chemistry box model MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere) was presented by Sander et al. (2005). The results of MOCCA (Sander and Crutzen, 1996) were reproduced by MECCA, and a comparison between these two box models was conducted. Furthermore, MECCA was used by Sander et al. (2006) to reveal the triggering of the bromine explosion and ozone destruction mechanism. It was proposed by Sander et al. (2006) that the precipitation of CaCO 3 from the freezing sea water is 5 a key factor of ODEs as it reduces the buffering ability of the sea water, which facilitates the acid-catalyzed activation process.
However, by utilizing the equilibrium thermodynamic model FREZCHEM, Morin et al. (2008) indicated the risk of using the results from Sander et al. (2006). Based on the finding of Dieckmann et al. (2008) that ikaite (CaCO 3 · 6 H 2 O) precipitates in brine instead of calcite (CaCO 3 ), Morin et al. (2008) suggested that if calcite is the major calcium carbonate mineral that precipitates, the results from Sander et al. (2006) are valid. On the contrary, if ikaite precipitates during the brine freezing 10 instead of calcite, the alkalinity of aerosol would not be buffered, which contradicts the conclusion made by Sander et al. (2006). Later on, in a joint publication of Sander and Morin (2010), the authors showed by re-analyzing the results from MECCA and FREZCHEM that the conclusions of Sander et al. (2006) hold as long as the mechanism of evapoconcentration during the formation of aerosol particles is considered. In addition, they also suggested that the conclusions of Sander et al. (2006) and Morin et al. (2008) can be reconciled by the introduction of a bromide/alkalinity ratio.

15
The one-dimensional model studies of ODEs started from Lehrer et al. (2004) by using FACSIMILE to inquire the relative importance of sea salt aerosols and the fresh sea ice surface for the ozone destruction. In their 1-D model, the vertical distribution of the turbulent diffusivity is parameterized using a linear function of height. It was indicated by their model that the primary bromine source during ODEs is possibly the fresh sea ice surface which provides a huge amount of reactive surface and nearly inexhaustible bromine compared to the sea salt aerosols. The prerequisites for the occurrence of ODEs were also 20 discussed by Lehrer et al. (2004). Another 1-D chemical transport model, THAMO (Tropospheric HAlogen chemistry MOdel), was developed by Saiz-Lopez et al. (2008) with the aim of clarifying the source strengths of bromine and iodine required for maintaining the vertical profiles of BrO and IO observed during the CHABLIS field campaign. In their model, prescribed Br 2 and atom I fluxes from the snowpack are added, and the heterogeneous reprocessing of bromine and iodine on sea salt aerosols is also included. It was found that an initial Br 2 flux is required for keeping the model simulations in accordance with 25 the DOAS measurements of the CHABLIS campaign. The recycling of bromine on sea salt aerosols was also found essential  (Bott et al., 1996;Bott, 1997Bott, , 2000. MISTRA is a boundary layer model initially developed for the investigation of the cloud-topped planetary boundary layer. Afterwards, the gas-phase and aqueous-phase reactions were implemented in MISTRA (von Glasow et al., 2002a, b;von Glasow and Crutzen, 2004). The new model is named as MISTRA-MPIC (MISTRA Max-Planck-Institut für Chemie) and used for the study of the marine boundary layer at mid-latitude. The halogen 35 chemistry in a marine boundary layer is modeled in MISTRA-MPIC, and the vertical structure of chemical species such as BrO in a cloud-free marine boundary layer is also captured. MISTRA was then modified by Piot and von Glasow (2008) to reproduce Arctic conditions for the investigation of the role of frost-flower derived aerosols, open leads and re-release processes on the snowpack in the destruction of ozone during ODEs. It was found that the recycling process on snow is the most important process for the existence of high-level bromine observed in the polar boundary layer. The effects brought about by the change 5 of meteorological parameters such as the temperature of the ambient air were also studied.
In another publication of Piot and von Glasow (2009) was found that significant bromide is still available in the liquid-like layer of the snowpack, which indicates that the snowpack is able to provide adequate reactive bromine to sustain the BrO level observed during the GSHOX experiment.
In order to address the influence of snowpack on the ozone loss in the polar boundary layer, Toyota et al. (2014) developed the 1-D model PHANTAS (PHotochemistry ANd Transport between Air and Snowpack). Therein, great attention was paid to the model configuration in order to realistically simulate the exchange between the snow interstitial air and the ambient 20 air in the boundary layer. Both the rapid ozone depletion in the boundary layer and the enhanced bromine in the interstitial air were successfully simulated applying PHANTAS. It was also found that in the top layers of the snowpack, the release of Br 2 is dominantly driven by the bromine explosion mechanism, whilst in deeper layers of snowpack, both the aqueous radical chemistry and the bromine explosion mechanism contribute substantially to the bromine evasion to the boundary layer.
The three-dimensional model studies on the tropospheric ODEs during polar spring started from Zeng et al. (2003). In their 25 3-D regional chemistry transport model (RCTM), fed with the data of BrO obtained from the spectral measurements aboard GOME, low ozone episodes occurring at northern high latitude during the TOPSE campaign were successfully captured.
Moreover, the correlation coefficients between the observed and modelled ozone temporal variations at different sites such as Alert and Barrow were found to be above 0. the heterogeneous reactions responsible for the reactivation of inactive bromine species in aerosols were implemented as well.
It was found that the bromine explosion cannot be predicted if the blowing snow module is switched off, which indicates the importance of blowing snow for the bromine explosion.
A major advancement of the 3-D models for studying bromine explosion and ozone consumption in the Arctic was made by Zhao et al. (2008). They used a global 3-D chemistry and transport model GEM-AQ/Arctic to investigate the spatial structure  (Sander et al., 2005(Sander et al., , 2006 described above was also implemented in GEM-AQ/Arctic, replacing the original chemistry module and solver. Moreover, the precipitation of CaCO 3 for triggering the bromine explosion outlined by Sander et al. (2006) is also included in their global model. A good agreement was achieved between the GEM-AQ/Arctic simulations and the satellite 10 observations. Additionally, it was shown in their model results that besides the halogen chemistry also the air temperature, atmospheric circulation, and long-range transport of pollutants make great contributions to the springtime ozone depletion in Arctic.
Another version of the 3-D model GEM-AQ with the incorporation of the gas-phase and heterogeneous bromine chemistry was employed by Toyota et al. (2011) to discover the surface sources of bromine for the occurrence of ODEs and BrO-clouds 15 observed in the boundary layer of the Arctic (Chance, 1998;Richter et al., 1998;Platt and Wagner, 1998;Wagner et al., 2001).
A source representing the oxidation of Brby ozone at the snow/ice-covered surfaces was included which later proved to be the primary source of reactive bromine. The simulation results compare reasonably well with ground-based measurements of the temporal behavior of surface ozone. The vertical distribution of ozone obtained by the model is also in consistence with the ozonesonde data measured at different locations (Alert, Resolute and Ny Ålesund). In 2013, Cao and Gutheil coupled the 3-D 20 model to a LES (Large Eddy Simulation) of the atmospheric flow using OpenFOAM in order to gain a better understanding of the turbulent mixing of ozone and bromine related species in the boundary layer during ODEs. In their paper they also discuss the dependence of the ozone depletion rate on meteorological conditions such as the wind speed and the boundary layer stability.
The above discussed theoretical studies revealed that ODEs are influenced by the joint effect of horizontal advection, vertical 25 convection, and local chemistry. The potential to provide the information which may enable to fully understand the underlying physicochemical processes is offered by the numerical simulation of ODEs of high dimensions. However, the possibility of conducting high-dimensional computations is usually limited by the effort to calculate the chemical production and destruction source terms in the governing equations. If a complex reaction mechanism is adopted in this kind of simulations, the execution time for the estimation of chemical source terms can exceed the one for solving transport equations by one order of magnitude. 30 Commonly, a one-step global reaction is adopted in 2-D or 3-D simulations of the atmospheric chemistry which leads to affordable computing times (Warnatz, 1992). However, the coefficients used in the one-step global reaction approach are usually empirical which cannot be estimated without the aid of experimental studies. An alternative treatment of the chemistry in numerical studies is to choose a couple of elementary reactions which are able to describe the important features of the reaction system. However, the selection of elementary reactions needs special attention. If the treatment of the chemistry in 35 the computation is too rough and arbitrary, the estimation of the chemical production and consumption would be inaccurate and, consequently, the simulation results may be far off from reality. Thus, an appropriate reaction mechanism with reduced size and adequate accuracy is required so that the multi-dimensional simulations of ODEs are applicable. The adopted reaction mechanism is required to maintain the major properties of the original complex reaction mechanism. Moreover, the size of the reaction mechanism should be small to enable the efficient computation of high-dimensional equations.

5
Previously, various approaches have been proposed to simplify a complex reaction mechanism while the major properties of the original reaction mechanism are maintained. The first type of the methods is represented by the rate spectrum analysis of the reaction scheme (Turányi et al., 1989). Therein, the rates of competing reactions are compared. However, this method is inapplicable for a very complex chemical system in which the influence of each reaction differs significantly on certain kinetic features of the system. The second type of methods is given by the grouping of the reactions into several categories 10 (Edelson and Allara, 1980). The reaction which makes negligible contribution to the flux within its category is considered to be unimportant, and, thus, can be eliminated from the reaction scheme. However, the grouping procedure of the reactions depends strongly on the experience of the investigator, which brings some arbitrariness to the application of this method. Moreover, the kinetic information inherent in the reaction mechanism may be lost if reactions belonging to a certain chemical kinetic structure are grouped separately. The third method is based on the analysis of the production rate of free radicals (Lifshitz and 15 Frenklach, 1975). Reactions which make small contributions to the production or destruction of free radicals are omitted. The fourth approach for obtaining the reduced mechanism is referred to as rate-of-production analysis or reaction path analysis (Glarborg et al., 1986). Herein, the individual contribution of each reaction to the overall production rate of selected chemical species is estimated. Also, the major formation pathway of the chemical species is identified.
A more complicated method which is able to extract the inherent information from the kinetic mechanism and discriminate 20 reactions from complex chemical reaction schemes is the sensitivity analysis (Rabitz et al., 1983;Turányi, 1990b). Sensitivity analysis is a routinely-used technique which reveals interdependence between the species concentrations and the change of reaction rates. Finally, it helps to reduce the size of the reaction mechanism. At present, the sensitivity analysis is often conducted in a "brute force" way (Dodge and Hecht, 1975;Valko and Vajda, 1984) by simply applying a perturbation of reaction rate coefficients by a fixed amount to see the corresponding deviations of particular species concentrations. This is time-consuming 25 and provides only a limited accuracy. Since the 1980s, more systematic and economical approaches for the estimation of the sensitivities have been proposed. The most straightforward one is the concentration sensitivity analysis (Turányi, 1990b), which is capable of displaying the relative importance of an individual rate coefficient on a group of species concentrations and identifying the rate-determining steps in the reaction scheme. The concentration sensitivity analysis has also been used for screening reactions from a complicated reaction mechanism in various fields such as the investigation of combustion (Dougherty and 30 Rabitz, 1980;Warnatz et al., 2001) and atmospheric chemistry (Turányi and Bérces, 1990).
However, the calculation of the concentration sensitivity coefficients depends on the time interval used in the model, which denotes that the concentration sensitivity is inherently time-depending. Thus, the accuracy of the sensitivity estimation is associated with the length of the time interval. Aside from this, the concentration sensitivity analysis only discovers the importance of an individual input parameter of the chemical kinetic system for several species concentrations while usually the concentra- 35 tions are influenced by a group of parameters. Thus, another mechanism reduction approach, namely the principal component analysis, is suggested (Vajda et al., 1985). The principal component analysis is performed based on the calculation of eigenvalues and eigenvectors of the matrix S T S, wherein S denotes an array of the relative concentration sensitivity matrices.
Strongly interacting reactions which significantly influence the concentration of particular species are identified by principal components with large eigenvalues. Moreover, further mechanistic details of the chemical kinetic system such as the species 5 in quasi-steady state can also be provided by the principal component analysis.
Although the mechanism reduction techniques mentioned above have been applied in the investigation of various atmospheric phenomena such as the chemistry in clouds (Pandis and Seinfeld, 1989), to date, little effort has been paid to the reduction of the reaction mechanism of ODEs using these techniques. Previous research of the author and co-workers (Cao and Gutheil, 2013;Cao et al., 2014) focused on this aspect to some extent. In these studies, relative concentration sensitivities of a 10 reaction mechanism of ODEs were computed, and the simplification of the mechanism based on the concentration sensitivities was prospected. However, the details of the reduction processes were not presented. Besides, the reduction technique discussed in these studies (Cao and Gutheil, 2013;Cao et al., 2014) is only limited to the concentration sensitivity analysis. Thus, in the present study, two different approaches, concentration sensitivity analysis and principal component analysis, were applied for reducing the size of a complex reaction mechanism of ozone depletion and halogen release. The results obtained by using 15 these two methods are compared and discussed. The connection between different chemical reactions is also revealed during the reduction processes.
The manuscript is organized as follows. In Sect. 2, the governing equations of the reduction approaches used in this study are presented. The criteria for screening unimportant reactions from the original reaction scheme are also given in this section.
Then the computational results of this work are shown in Sect. 3. Reduced reaction mechanisms derived after the analyses 20 are presented and compared as well. The species in quasi-steady state is also revealed by the implementation of the principal component analysis. At last, major conclusions made in the present study are addressed in Sect. 4. Further simplification of the reduced reaction mechanism is also prospected.

Mathematical Model and Methods
A complex chemical reaction system can be denoted as with the initial condition c| t=0 = c 0 , where c is a column vector of species concentrations. k in Eq.
(1) is a vector of reaction rate coefficients and t denotes time. E represents the source term of local surface emissions. In the present study, a reaction mechanism with the involvement of bromine containing compounds and nitrogen related species is adopted from the previous box model study (Cao et al., 2014). Furthermore, the reaction rate constants are updated with the latest chemical kinetic data 30 (Atkinson et al., 2006). Several photolysis reactions are also added to complete the previous reaction mechanism. The reaction scheme used in the present study consists of 39 chemical species and 92 reactions and is listed in the supplementary material.
In this chemical reaction mechanism, bromine is considered as the only halogen species while the chlorine related species and reactions are excluded. The influences brought about by the inclusion of chlorine related species and reactions have been discussed by Cao et al. (2014). As the focus of the present study is laid on the reduction of the reaction mechanism instead of pursuing a more precise prediction of the temporal change of each chemical species, the choice of this chlorine-free reaction mechanism does not affect the metrics of the present study.

5
As mentioned above, two different approaches, concentration sensitivity analysis and principal component analysis, are followed to reduce the size of the reaction mechanism of ODEs. The governing equations of these two approaches and the procedures of the mechanism simplification are presented below.

Concentration Sensitivity Analysis
The importance of the j-th reaction for the i-th chemical species is depicted by the relative concentration sensitivity S ij , which 10 can be written in the form of In Eq. (2), c i is the concentration of the i-th chemical species, and k j denotes the rate coefficient of the j-th reaction. S ij = ∂c i /∂k j is the absolute concentration sensitivity, and the unit of S ij depends on the order of the j-th reaction. In order to compare the sensitivity coefficients belonging to different reactions, the normalized version of the sensitivity coefficient, S ij , 15 is introduced by multiplying S ij with k j /c i . The obtained relative concentration sensitivity S ij is thus a dimensionless variable which represents the percentage of change in the i-th species concentration due to a 1% change of the j-th reaction rate constant. The evaluation of the relative concentration sensitivity is helpful for discovering the interdependence between the solution of Eq. (1) and the input parameters such as the reaction rate constants k j .
The calculation of S ij is performed by differentiating the i-th component of Eq.
(1) with respect to k j . With the assumption 20 that the local emissions are independent of the rate constants, Eq. (1) becomes The upper limit of the sum in Eq. (3), ns, denotes the total number of the chemical species included in the model. By substituting ∂c i /∂k j and ∂c l /∂k j in Eq.
(3) with the absolute concentration sensitivities S ij and S lj , we obtain another form of Eq. (3) as follows, The second term of the right hand side of Eq. (4) denotes the direct influence on the concentration of i-th species caused by the change of the j-th reaction rate constant. Moreover, as a result of this direct variation in the j-th reaction rate, indirect effects on the concentrations of other species are induced via the coupled kinetic system, which contributes to the solution of Eq. (4).
This indirect effect of parameter change is indicated by the first term of the right hand side of Eq. (4). After obtaining the solution S ij of Eq. (4), the relative concentration sensitivity can be calculated by multiplying S ij with k j /c i . The computation of the absolute concentration sensitivity is conducted by using the subroutine SENS in the chemical kinetic software KINAL (Turányi, 1990a). The Decomposed Direct Method (Valko and Vajda, 1984) is implemented in KINAL for solving Eq. (4), which has been proved robust and efficient (Turányi, 1990b).
The concentration sensitivity analysis is a useful measure of how sensitive a specified species concentration is to a par-5 ticular reaction rate constant. Thus, reactions with large absolute values of sensitivities are identified as important and ratedetermining. For the purpose of deriving an appropriate reduced reaction mechanism of ODEs, it is needed to remove the least important reactions from the original reaction scheme. In the present study, we consider the j-th reaction as unimportant if the 10 is fulfilled. In Eq. (5), ns is the total number of the chemical species as mentioned above, and nt denotes the total time step in the computation. The criterion shown in Eq. (5) means that the j-th reaction is considered as unimportant if its relative concentration sensitivity for all species at any time point is smaller than 10% and then can be eliminated from the system.
Although the concentration sensitivity analysis is powerful for constructing a minimal chemical reaction set describing ODEs, it has some disadvantages. As the concentration sensitivity represents the response of a system to a perturbation of the 15 reaction rate at an earlier time point, the concentration sensitivity coefficient is inherently time-depending, and its estimation depends on the adopted time interval. Thus, the concentration sensitivity obtained by solving Eq. (4) corresponds to a time interval instead of a fixed time point. A drawback brought about by this inherent time-depending property of the concentration sensitivity is the "memory effect" (Vajda and Turányi, 1986;Turányi et al., 1989). Within the time interval for the integration of sensitivity coefficients, if a reaction is indicated as important at any stage of the time interval, it is identified as important during 20 the whole time interval since this property is maintained in the concentration sensitivity analysis. As a result, this reaction cannot be eliminated from the original reaction scheme. In addition, although the concentration sensitivity analysis is capable of clarifying the dependence of species concentrations upon a particular rate coefficient, usually the species concentrations are influenced by a combination of input parameters such as a ratio of two reaction rate constants. This interdependence between the species concentrations and the groups of parameters can hardly be revealed by the concentration sensitivity analysis. 25 Therefore, in the present study, another mechanism reduction approach is chosen, namely the principal component analysis, and applied to the original reaction mechanism of ODEs, which is discussed in the section below.

Principal Component Analysis
The principal component analysis provides an effective means of screening unimportant reactions from a complex reaction scheme so that a tractable reaction mechanism can be derived. To perform the principal component analysis, we first introduce 30 the normalized rate parameter, α, of which the j-th component can be expressed as α j = ln k j , j = 1, ..., np.
Herein, np denotes the total number of the chemical reactions considered in the model.
The definition given by Eq. (6) enables to express the response of the reaction mechanism due to a variation in the rate parameters, α, using the response function Q(α) as follows, In Eq. (7), the subscript "i, m" represents the concentration of i-th species at the m-th time point. ns and nt are the total 5 number of the chemical species and the time step, respectively. α 0 denotes a matrix of the normalized parameters with original values.
By applying a Taylor series expansion (Vajda et al., 1985) and a Gauss approximation (Bard, 1974) on Eq. (7), the response function Q ( α) is approximated as 10 Herein, Q(α) is an approximate response function. ∆α denotes a column vector of the parameter variations which has the j-th component, S in Eq. (8) represents an array of the relative concentration sensitivities corresponding to the time instance t(m), m = 1, . . . , nt, and the m-th element of S has the form of S is, thus, a matrix with a dimension of (nt · ns) × np. The eigenvalue-eigenvector decomposition technique is then applied on the matrix S T S in Eq. (8), which leads to where U is a matrix of normed eigenvectors of S T S. Λ represents a diagonal matrix consisting of the eigenvalues λ 1 , λ 2 , ..., λ np 20 of S T S. By replacing S T S in Eq. (8) with the expression shown in Eq. (11), we obtain the approximate response function in the form of We then define the principal component Ψ = U T α in order to write U T (∆α) in the above equation as a variation of the principal component, As a result of the definition of the principal component, Eq. (12) becomes wherein ∆Ψ j 2 = (∆Ψ j ) T (∆Ψ j ). From Eq. (14), it can be seen that the eigenvalue λ j depicts the significance of a group of reactions for the change of the system. ∆Ψ j in Eq. (14)  Therefore, after obtaining the relative sensitivity coefficient S ij of the reaction mechanism representing the chemistry of 10 ODEs, principal component analysis is performed on the matrix S T S. In that way np eigenvalues and the associated eigenvectors are obtained. In order to derive a reaction mechanism with reduced size, we remove the reactions belonging to the principal components with small eigenvalues. It has been proved in previous studies (Vajda et al., 1985) that if the eigenvalue λ j is smaller than ns × nt × 10 −4 , the variation in the concentration of each chemical species in the system is less than 10% at each instance in time. This selection criterion is also adopted in the present study. Aside from this, in the principal components 15 corresponding to large eigenvalues, elements with values less than 0.2 are also considered as unimportant and, thus, can be removed from the original mechanism. It was estimated that these elements contribute less than 4% to the total variation of concentrations (Vajda et al., 1985).
It should also be noted that by removing the j-th reaction from the original chemical system, we actually define the reaction rate constant of the j-th reaction k j = 0, for which Eq. (6) is invalid. To solve this problem the parameter is introduced (Vajda et al., 1985). It can be easily deduced that α j and the parameter α j discussed above yield the same normalized sensitivity S ij . Therefore, the investigation of S T S for α j described in the present manuscript is also valid for α j in the simplification process of the reaction mechanism.
Thus, at first, the original reaction mechanism described above is implemented in the box model KINAL (Turányi, 1990a) to 25 capture the temporal behavior of ozone and principal bromine species. The boundary layer height used in the model is defined as 200 m. This is considered to be a representative value since the boundary layer height under typical polar conditions ranges from 100 m to 500 m (Stull, 1988). The initial atmospheric composition used in the model is taken from the previous box model study (Cao et al., 2014), which is listed in Tab. 1. In the present model, the prescribed 0.3 ppt Br 2 and 0.01 ppt HBr play the role of triggering the bromine explosion mechanism and the associated ozone consumption. Emissions of nitrogen 30 oxides (NO x ), HONO, H 2 O 2 and HCHO from the underlying surface are also implemented in the model according to the measurements (Jones et al., 2000(Jones et al., , 2001Jacobi et al., 2002;Grannas et al., 2007), and listed in Tab. 2. The ratio of HONO and NO 2 is assumed to be unity (Grannas et al., 2007).
The photolysis reaction rates are estimated by using a three-stream radiation transfer model (Röth, 1992(Röth, , 2002 with an assumption of SZA = 80 • and a surface albedo unity. At present, the daily variation in SZA and the dark reactions are not accounted for in the model. The influences brought about by the inclusion of the change in SZA and the dark reactions have been proved to be small in previous investigations (Lehrer et al., 2004;Cao et al., 2014). The heterogeneous reactions representing the bromine recycling processes on the surfaces of the ice/snowpack and the suspended aerosols are also included 5 in the original reaction mechanism of ODEs. It has been proved that the rates of these heterogeneous reactions critically control the bromine amount in the boundary layer and the depletion rate of ozone during ODEs (Cao et al., 2014). In the present study, the parameterizations of these heterogeneous reactions are also adopted from the box model study by Cao et al. (2014). Recently, a snowpack module representing the mass exchange between the snowpack and the ambient air has been developed by the author and co-workers (Cao et al., 2016b). However, this snowpack module was not applied in the reaction 10 mechanism used in the present study.
Later, the concentration sensitivity analysis is applied on the original reaction mechanism of ODEs to identify the most influential reactions during different time periods. The temporal behavior of the concentration sensitivity coefficient for each reaction is also captured. Then the reactions with a maximum absolute value of the relative concentration sensitivity less than 10% are removed from the original reaction scheme so that a reduced reaction mechanism is obtained. 15 After the computation of the relative concentration sensitivity coefficients, the principal component analysis is performed on the matrix S T S. Principal components corresponding to small eigenvalues are regarded as making negligible contributions to the overall response of the whole system and, thus, can be eliminated. Since 38 chemical species (all species in the mechanism excluding N 2 ) and 94 time steps are focused on in the present simulation, the criterion of eigenvalues used for dividing important and unimportant principal components is calculated as 38 × 94 × 10 −4 = 0.3572. Moreover, as discussed above, if 20 an element in a relatively important principal component has a value of less than 0.2, this element can be also removed due to its relatively small contribution to the total variation of the system.
In the following section, the most important computational results are presented and discussed.

Results and Discussion
At the beginning of this study, we ran the box model KINAL with the implementation of the original complex reaction scheme 25 to capture the temporal change of principal chemical species within the 200 m boundary layer. Figure 1 displays the development of the mixing ratios of ozone and principal bromine species with time. According to previous studies (Cao et al., 2014), the whole ODE can be divided into three periods. In the first period named "induction stage" which corresponds to the time period before day 3 in Fig. 1, the consumption rate of ozone in the boundary layer is less than 0.1 ppb h −1 . Moreover, the mixing ratio of ozone remains steady at a background level (∼ 40 ppb). Due to the presence of ozone in this time period, the 30 formed Br atoms are rapidly oxidized to be BrO, and thus can be hardly observed at this time stage. In contrast to that, the mixing ratios of HOBr and BrO steadily increase within this time period, which makes them the major bromine containing compounds.
After day 3 (see Fig. 1), as bromide is continuously activated from the ice/snow surface via the bromine explosion mechanism, the total bromine loading in the ambient air keeps increasing. As a result, the depletion rate of ozone exceeds 0.1 ppb h −1 , and the mixing ratio of ozone declines rapidly until a value lower than 10% of its original level (4 ppb) is achieved, on approximately day 4.6. This period (from day 3 to day 4.6) is thus named "depletion stage". Within the depletion stage, the mixing ratios of BrO and HOBr reach their peaks and then drop instantaneously to a level less than 1 ppt due to the nearly complete 5 removal of ozone in the boundary layer.
After the depletion stage, the last time period of the ODE is the "end stage" in which the ozone mixing ratio is lower than 4 ppb. At the beginning of the end stage, Br atom becomes the most dominant brominated species and its mixing ratio increases to a peak value of approximately 170 ppt (see the black solid line in Fig. 1). This large amount of Br is then absorbed by the aldehydes in the atmosphere, and the total bromine amount left in the ambient air becomes steady towards the end of the 10 simulation (see the purple solid line in Fig. 1). After day 6, the dominant bromine species in the ambient air is mostly in the form of HBr, which is in consistence with the previous finding (Langendörfer et al., 1999). The features of the temporal change of these chemical species are similar to those obtained in previous studies (Lehrer et al., 2004;Cao et al., 2014Cao et al., , 2016a) except a slight difference in the onsets of the aforementioned time periods. The discrepancy between the results from different model studies is attributed to the extension of the previous reaction mechanism in the present study. 15 After capturing the temporal behavior of principal chemical species, the original reaction mechanism consisting of 39 species and 92 reactions is processed with the concentration sensitivity analysis, which is shown in the next subsection.

Concentration Sensitivity Analysis of the Reaction Mechanism of the ODE
The relative concentration sensitivity coefficient of each chemical species for all the reactions in the original mechanism within each time step is calculated by performing the concentration sensitivity analysis. Figure 2 depicts the sensitivity coefficients 20 of ozone and BrO for all the reactions in the mechanism within the time interval [day 3.9, day 4] which resides within the depletion stage. It can be seen that the sensitivity of the mixing ratios of ozone and BrO to each chemical reaction in the mechanism is clearly indicated by the values of the sensitivity coefficients.
In Fig. 2(a), reaction (R15), which denotes the heterogeneous bromine activation from the ice/snow-covered surfaces, is identified as the most important reaction for both ozone and BrO mixing ratios at this time stage. This is related to the fact that 25 during the depletion stage, the HOBr concentration increases vigorously until its peak value is achieved. Thus, a large amount of bromide is activated from the ice/snow surfaces. As a result, reaction (R15) critically determines the total bromine amount in the air and the depletion rate of ozone in this time period. Due to the same reason, reactions with the involvement of HOBr, (R10) and (R11), play significant roles in influencing the mixing ratios of ozone and BrO, which is indicated in Fig. 2(a).
Aside from the reactions associated with HOBr, reactions which are able to produce or consume Br atoms are also influential since Br reacts with ozone directly and rapidly. Thus, the significance of the Br associated reactions, (R5), (R7), (R17) and (R20), is displayed in Fig. 2(a) according to the relatively larger absolute values of the sensitivity coefficients. By evaluating the sensitivities at an earlier time instance (not shown here), we found that the Br associated reactions are the most dominant reactions for ozone and BrO in the induction stage. However, from Fig. 2(a), it can be concluded that during the depletion stage, the importance of the HOBr associated reactions exceed the ones with the involvement of Br atoms.
In Fig. 2(b), it is shown that the nitrogen related reactions are less important for the mixing ratios of ozone and BrO in the boundary layer compared to the bromine related reactions. However, the moderate role of the hydrolysis reaction of BrONO 2 , (R90), is depicted. The reason is attributable to the formation of HOBr in this reaction, which may also speed up 5 the bromine activation during this time period. Moreover, we found that the reactions with the involvement of BrONO 2 are of more importance during the end stage, which is consistent with the previous box model study (Cao et al., 2014).
The results discussed above show that the values of the concentration sensitivity coefficients vary with time. This reveals that the importance of each reaction differs significantly during different time periods. In order to clarify the temporal behavior of the concentration sensitivities, we also captured the development of the sensitivities of ozone for the dominant reactions with 10 time (see Fig. 3). As similar as the results shown in Fig. 2, it is seen in Fig. 3 that the HOBr associated reactions, (R10), (R11) and ( ODEs. The reduced reaction mechanism obtained after the simplification process, thus, contains 39 species and 81 reactions in total.
We applied the reduced reaction mechanism in the box model KINAL to see the deviations from the original reaction scheme. Figure 4 shows a comparison of the mixing ratios of ozone and principal bromine species computed by using the original reaction mechanism and the reduced reaction mechanism. It is found that both simulations give nearly identical results. The 25 depletion stage of the ODE predicted by the reduced reaction mechanism lasts from day 3 to day 4.6 which is consistent with the results obtained by using the original reaction scheme. The maximum deviations of the mixing ratios of Br, HBr, HOBr, BrO and Br tot between the original reaction mechanism and the reduced reaction mechanism are 0.6%, 0.5%, 0.6%, 0.1% and 0.2%, respectively. Thus, it can be seen that the calculated concentration changes of ozone and principal bromine species from these two mechanisms agree very well, within deviations of 1%.

30
The principal component analysis was then performed on the sensitivity matrix S T S of the original reaction scheme to simplify the mechanism, which is shown in the next subsection.

Principal Component Analysis of the Reaction Mechanism of the ODE
After obtaining the concentration sensitivities of the original reaction mechanism of the ODE, the matrix S T S is constructed and applied in the principal component analysis. The eigenvalues and the corresponding principal components of S T S are thus obtained which are listed in Tab. 4. The selection criterion λ i < 0.3572 is adopted for removing reactions from the scheme which has been discussed above. Moreover, in the principal components with large eigenvalues, if an element contributes less

10
It is found that the reactions identified as unimportant by using the sensitivity analysis are all covered by the principal component analysis. Moreover, 9 extra reactions to be eliminated from the reaction scheme are also revealed by the principal component analysis. The reason for the difference between these two approaches is possibly that in a complex reaction mechanism, usually the removal of one single reaction would cause significant variations in the temporal change of many chemical species while the elimination of a group of reactions may have only minor impact on the response of the system. In the con- 15 centration sensitivity analysis, the association between a particular species concentration and each reaction rate is clarified, which is suitable for screening reactions separately. In contrast to that, the principal component analysis is able to identify the dependence of the species concentrations on a group of reactions. Therefore, by performing the principal component analysis, we were able to remove more reactions from the original reaction mechanism. At last, the reduced reaction mechanism after the implementation of the principal component analysis consists of 39 species and 72 reactions (see Tab. 3).

20
After performing the principal component analysis, the reduced reaction mechanism was introduced into the box model to simulate the mixing ratio change of ozone and principal bromine species (see Fig. 5). It is seen that the results obtained by using the reduced reaction mechanism agree well with those simulated by using the original reaction mechanism. The depletion stage in simulation results using the reduced reaction mechanism starts on day 2.9 and finishes on day 4.4, which occurs a little bit earlier than that using the original reaction mechanism. The maximum deviations of Br, HBr, HOBr, BrO and Br tot are 6.6%, 25 1.6%, 8.9%, 2.5% and 5.2%, respectively. Thus, the variations in the mixing ratios of the principal bromine species considered in the model caused by the removal of 20 redundant reactions are less than 10% which is restrained by the selection criterion adopted. Therefore, the reduced reaction mechanism can satisfactorily capture the temporal change of each chemical species so that the requirement of the accuracy of the reaction mechanism is fulfilled for multi-dimensional computations of ODEs.
The principal component analysis was applied to different time periods of ODEs (induction stage, depletion stage and end rates of (R1) and (R3), which can be expressed as Herein, k 1 and k 3 are the rate constants of reactions (R1) and (R3), respectively. In order to confirm this finding, we multiplied the rate coefficient of both (R1) and (R3) by 100, the deviations of the system after this change are shown in Tab. 5. It is found that the peak values of the principal bromine species calculated from the reduced reaction mechanism and from the reaction 25 scheme with the modified reaction rates agree well within 1 percent. In contrast to that, if we only increase the rate constant of (R1) to 100 times of its original value while the rate constant of (R3) remains the same, the difference between the results from the original mechanism and the modified mechanism becomes larger (see Tab. 5). The maximum deviation amounts to 172% for the species HOBr. This finding also confirms the validity of the principal component analysis.
As O 1 (D) was found in quasi-steady state during ODEs, it is possible to directly estimate the mixing ratio of O 1 (D) at each 30 time instance according to the mixing ratios of other chemical species. It can be easily deduced that the mixing ratio of O 1 (D) during the depletion event is expressed as in which [O 3 ] and [N 2 ] are the mixing ratios of ozone and nitrogen, respectively. Thus, in multi-dimensional model studies of ODEs, aside from using a reduced reaction mechanism to shorten the time cost for the calculation of the chemical source terms, the species equation corresponding to O 1 (D) can be also neglected, which improves the efficiency of the computations.
However, the simplification of the reaction mechanism by using the quasi-steady state approximation is beyond the scope of the present study.

4 Conclusions and Future Developments
In the present study, two reduction approaches, namely the concentration sensitivity analysis and the principal component analysis, were applied on a reaction mechanism representing the chemistry of ODEs. The former was performed based on the ratio of the relative change of species concentrations and the variation of a particular reaction rate. The significance of each reaction in the original mechanism for various chemical species was identified. It was found that during the depletion of ozone, 10 reactions associated with HOBr are the most influential reactions as they determine the total bromine loading in the atmosphere within this time period. By removing 11 reactions which exhibit peak absolute values of the sensitivity coefficients lower than 10%, a reduced reaction mechanism of ODEs was derived. The difference between the results obtained by using the original reaction mechanism and the reduced reaction mechanism is negligible, which validates the sensitivity analysis.
The principal component analysis, on the other hand, is able to extract inherent information from a chemical kinetic system 15 via the eigenvalue decomposition of the matrix S T S. In the herein presented study, the reactions which have principal components with eigenvalues lower than a critical criterion were eliminated from the original reaction mechanism. Additionally, the elements which occupy less than 0.2 of each principal component were also removed. In total, 20 reactions were identified redundant using the principal component analysis and, thus, screened out from the original ODE mechanism. As a result, a reaction scheme with a reduced size consisting of 39 species and 72 reactions was obtained. The deviations of the mixing ratios 20 of ozone and principal bromine species between the original reaction mechanism and the reduced reaction mechanism after the principal component analysis are within 10%. This proves the suitability of the obtained reduced reaction mechanism for multi-dimensional simulations of ODEs. Apart from this, displayed by the principal components with the smallest eigenvalues, the chemical species O 1 (D) is identified to be in quasi-steady state. This facilitates a further improvement of the numerical efficiency in high-dimensional simulations of ODEs. 25 It is obvious that the applications of the concentration sensitivity analysis and the principal component analysis are not restricted to the reaction mechanism of ODEs which is the scope of the present study. Any chemical kinetic system which is expressible through a form of Eq. (1) can be analyzed and optimized by using these two approaches. Furthermore, the readers might have noticed that the two reduction approaches used in the present study do not lead to the removal of redundant species from the detailed reaction mechanism. The reason is that they, rather than employing erroneous simplifications (Vajda 30 et al., 1985;Turányi, 1990b), focus on all occurring chemical species. This, however, prohibits the identification of redundant species. Therefore, it seems worth to utilize another approach in order to further reduce the reaction mechanism of the ODE after the analyses performed herein. At present, three types of reduction methods are normally used in previous studies (Turányi and Tomlin, 2014): species lumping, timescale-separation method, and tabulation approaches. As discussed above, the quasisteady state approximation, which belongs to the category of the timescale-separation method, is capable of determining the concentrations of the species in quasi-steady state through a function of other species concentrations, which helps to remove the species in quasi-steady state from the kinetic reaction scheme. Apart from this, currently, another type of a timescale-based analysis approach, the so-called intrinsic low-dimensional manifolds (ILDM) (Maas and Pope, 1992), is being implemented by      The solid curves denote the simulation results obtained by using the original reaction mechanism, and the dashed curves represent the results for the reduced reaction mechanism which is obtained after the concentration sensitivity analysis. The curves of ozone mixing ratios use the axis of time while the principal bromine species correspond to the ozone coordinate axis. Note that the axis of time is in a reverse direction for clarity. The solid curves denote the simulation results obtained by using the original reaction mechanism, and the dashed curves represent the results for the reduced reaction mechanism which is obtained after the principal component analysis. The configuration of this figure is similar to that used in Fig. 4.