J Galaz
email: jdgalazm@gmail.com

R Cienfuegos

A Echeverria

S Pereira

C Bertin
email: celbertin@gmail.com

G Prato

J C Karich

Integrating tsunami simulations in web applications using BROWNI, an open source clientside GPU-powered tsunami simulation library

Keywords: Tsunami Simulation library Javascript GPU Web Visualization Conceptualization, methodology, software (BrowNI and TsunamiLab), validation, formal analysis, investigation, data curation, visualization, project administration. R.C.: Conceptualization, methodology, validation, resources, supervision. A Conceptualization, methodology, supervision Conceptualization, software, TsunamiLab

BROWNI is a WebGL -Javascript tsunami simulation

Introduction

Tsunamis have become one of the deadliest natural hazards in the world, causing more than 200,000 casualties and billions of dollars in economic losses in the last 20 years [START_REF] Kânoğlu | Tsunamis: bridging science, engineering and society[END_REF], which largely justifies the important efforts in increasing tsunami preparedness and awareness that different governments have undertaken over the last decades. Tsunami early warning systems, such as those from Japan, Chile, and the USA [START_REF] Catalán | Sistema integrado de predicción y alerta de tsunamis (sipat)[END_REF][START_REF] Titov | Development of MOST for real-time tsunami forecasting[END_REF][START_REF] Kamigaichi | Tsunami forecasting and warning[END_REF], are examples of these efforts which have enabled communities to react more quickly and in turn prevent casualties. These systems are based on state of the art scientific knowledge and include numerical modeling as an essential tool for performing wave propagation forecast [START_REF] Kânoğlu | Tsunamis: bridging science, engineering and society[END_REF].

Thanks to the recent advances in computational modeling, it is now possible to simulate complex scenarios with increased resolution (e. g., GeoClaw [START_REF]Clawpack software[END_REF], [START_REF] Berger | The GeoClaw software for depth-averaged flows with adaptive refinement[END_REF], COMCOT [START_REF] Wang | User manual for comcot version 1.7[END_REF], CoulWave [START_REF] Lynett | Modeling wave generation, evolution, and interaction with depthintegrated, dispersive wave equations coulwave code manual[END_REF], NeoWave [START_REF] Yamazaki | Proceedings and results of the 2011 NTHMP model benchmarking workshop[END_REF], EasyWave [START_REF] Christgau | A comparison of cuda and openacc: accelerating the tsunami simulation easywave[END_REF], ANUGA [START_REF] Nielsen | Hydrodynamic modelling of coastal inundation[END_REF], and MOST [START_REF] Titov | Development of MOST for real-time tsunami forecasting[END_REF]). However, gathering and analyzing data iteratively using numerical modeling is still costly owing to expensive hardware requirements, complex installation processes and heavy tasks of data pre-and post processing.

The latter has so far limited the development of innovative applications for disaster risk management that could integrate numerical simulations with geographical information systems (GIS) and visualization tools [START_REF] Merati | Integration of tsunami analysis tools into a gis workspace-research, modeling, and hazard mitigation efforts within noaa's center for tsunami research[END_REF]. Indeed, new visualization and interactivity options have not been fully exploited yet, even if some interesting efforts to commuicate tsunami risks to emergency managers and broad audiences have been made taking advantage of such capabilities. For example, Keon et al. [START_REF] Keon | Web-based spatiotemporal simulation modeling and visualization of tsunami inundation and potential human response[END_REF] integrated inundation estimation with agent-based modeling to develop a framework that allows for interactive visualization and exploration of human response. Similarly, Hsieh et al. [START_REF] Hsieh | Parallel tsunami simulation and visualization on tiled display wall using opengl shading language[END_REF] showed an application that allows for efficient computation and visualization of tsunami simulations on tiled displays, taking advantage of the synchronized computing power of a GPU cluster which reaches high performance while efficiently rendering its results. Most remarkably, the Center for Tsunami Research (NCTR) of the National Oceanic and Atmospheric Administration (NOAA) [START_REF] Titov | Development of MOST for real-time tsunami forecasting[END_REF] fosters international collaborations providing emergency managers with different tools to help assessing the hazard impact visualizing demographics, evacuation routes, infrastructure, and other variables that are relevant for the decision.

These applications can be related to a generalized need for new geovisualizations that can help in communicating hazards and risks to non-expert audiences, such as scientists from other fields, emergency managers, coastal communities, policy makers [START_REF] Teeuw | Free or low-cost geoinformatics for disaster management: Uses and availability issues[END_REF], as well as children, teachers, and adults in general [START_REF] Kinzel | Using educational tools and integrative experiences via geovisualizations that incorporate spatial thinking, real world science and ocean literacy standards in the classroom: a case study examined[END_REF]. For example, regarding flood hazards, Jacquinod et al. [START_REF] Jacquinod | Automated production of interactive 3d temporal geovisualizations so as to enhance flood risk awareness[END_REF] created three-dimensional (3D) visualizations using GIS to increase awareness, and [START_REF] Curebal | Gis-based approach for flood analysis: case study of keçidere flash flood event (turkey)[END_REF] showcased how GIS can be useful for analyzing inundation areas by integrating it with engineering software tools such as HEC RAS [START_REF] Brunner | HEC-RAS river analysis system: hydraulic reference manual[END_REF].

GPU computing, on the other side, has allowed users to accelerate computations at a lower cost when compared to traditional central processing unit (CPU) clusters; GPU also produces complex visualizations efficiently. Recent examples of these are the open source software packages TsuPy [START_REF] Schäfer | Tsupy: computational robustness in tsunami hazard modelling[END_REF] and Celeris [START_REF] Tavakkol | Celeris: A gpu-accelerated open source software with a boussinesq-type wave solver for real-time interactive simulation and visualization[END_REF]. TsuPy is a script-based Python library that uses the NVIDIA Compute Unified Device Architecture CUDA [START_REF] Nickolls | Scalable parallel programming with cuda[END_REF] to calculate wave propagation and inundation at the regional scale with a shallow water solver, obtaining high-performance computations on NVIDIA GPUs. Celeris is a Boussinesq wave model that uses Microsoft Direct3D API [START_REF] Blythe | The direct3d 10 system[END_REF] shaders to calculate and render the simulations in a 3D interactive view with a graphical user interface that allows users to change parameters on the fly. As reported in [START_REF] Tavakkol | Celeris base: An interactive and immersive boussinesq-type nearshore wave simulation software[END_REF] Celeris was recently ported to the game engine Unity3D to go even further into producing an inmersive virtual reality user interface. These tools demonstrate that it is possible to bridge the gap between modellers (expert users) and decision makers by facilitating the tasks of software and hardware setup, configuration of scenarios, data preparation, post-processing, visualization and comparison with other sources of evidence.

The TsuPy software requires specific vendor hardware for its execution, complex installation processes and dependency management, and the usage of server-side computations in order to be used from an interactive interface; and while Celeris can take advantage of Unity3D being multiplatform, adding new interactions and sources of data requires editing its source code. However, another approach not yet explored in the literature is to take advantage of web-browsers to run the computations. Among other things, applications running on web browsers benefit from an almost null cost of installation and plenty of tools for developers to create graphical interfaces where users could gather and analyze data from simulations and other sources more efficiently. Moreover, browsers can also benefit from GPU computing to speedup calculations using the WebGL application programming interface (API) [START_REF] Marrin | Webgl specification[END_REF]. WebGL is a web standard with a JavaScript API designed for efficiently rendering 2D and 3D interactive graphics, and currently, the only interface that can do this without any plugins or browser extensions. This gives web browsers the advantage of directly mixing simulation results with other user-interface elements, handling asynchronous user interactions and data updates.

In this study, we examine the capabilities of GPU computing from the web browser to reduce the cost of simulations and take advantage of interactivity and visualization to improve tsunami risk communication and scientific outreach.

We introduce a tsunami-simulation library, the Browser's Numerical Interface (BROWNI), which has an API designed to facilitate the creation of interactive tsunami simulations and visualizations directly on the web browser. In section 2, we give an overview of the mathematical equations and numerical algorithm, to give the relevant context to understand the implementation of the library. Then we explain the software implementation of BROWNI in section 3, where we also discuss some limitations of this approach. In section 4 we validate the results of BROWNI by comparing the timeseries of wave heights against two real scenarios (Tohoku, 2011 and Chile, 2010), demonstrating its capabilities to obtain accurate results, while making efficient use of the available graphics card (dedicated or integrated). Then in section 5 we give two examples of scientific-outreach software that were built with BROWNI: TsunamiLab, a web application available from www.tsunamilab.cl, and the TsunamiLab-Pool, an interactive augmented-reality version of TsunamiLab designed for public exhibitions.

Mathematical model

The mathematical modeling of tsunami waves is a multi-scale problem that usually requires coupling several algorithms to accurately cover all the physical phenomena observed from wave generation and propagation in the ocean, to wave shoaling and run-up [START_REF] Leveque | Tsunami modelling with adaptively refined finite volume methods[END_REF][START_REF] Wang | User manual for comcot version 1.7[END_REF]27]. Since the focus of our work is on interactive applications for tsunami-risk management and communication, we implemented in BROWNI a far-field tsunami model that has been already used by others in operational settings [27,[START_REF] Christgau | A comparison of cuda and openacc: accelerating the tsunami simulation easywave[END_REF][START_REF] Titov | Development of MOST for real-time tsunami forecasting[END_REF][START_REF] Titov | Development of MOST for real-time tsunami forecasting[END_REF][START_REF] Imamura | Tsunami modelling manual (ver-3.1.4)[END_REF]. The implementation of more complex regimes and algorithms is left for future work.In this section we only recall these equations and its numerical resolution without further developments or improvements.

Model equations

In the ocean, tsunami waves have a characteristic wavelength ≈ 100 km, a characteristic amplitude ≈ 1 m, and a characteristic water depth of ℎ 0 ≈ 4 km, which means that waves are relatively long (ℎ 0 ∕ ≪ 1) and of small amplitude ( ∕ℎ 0 ≪ 1) [START_REF] Dean | Water wave mechanics for engineers and scientists[END_REF]. Under these characteristics, as presented in [START_REF] Liu | Numerical simulations of the 1960 chilean tsunami propagation and inundation at hilo, hawaii[END_REF], a suitable approximation is given by the linear shallow water equations in spherical coordinates:

+ 1 ( ) + ( cos ) = 0 + ℎ cos = + ℎ = - (1)
where is the time coordinate; , and are the longitude and latitude geographical coordinates, respectively; , , and are the wave height and longitudinal and latitudinal momentum components, respectively; = 6378 km is the earth's radius; = 9.81 m/s 2 is the earth's gravitational acceleration; ℎ( , ) is the bathymetry at location ( , ), where ℎ > 0 indicates underwater floor; and = 2 sin( ) is the Coriolis factor with = 7.29 × 10 -5 [ ∕ ], which is the earth's rotation frequency.

Numerical discretization

Similar to [START_REF] Christgau | A comparison of cuda and openacc: accelerating the tsunami simulation easywave[END_REF], [START_REF] Wang | User manual for comcot version 1.7[END_REF], and [27], equation ( 1) is discretized using finite differences with a second order in space and time leapfrog scheme as follows:

+1∕2 , - -1∕2 , Δ + 1 cos +1∕2, --1∕2, Δ + , +1∕2 cos +1∕2 -, -1∕2 cos -1∕2 Δ = 0 +1 +1∕2, -+1∕2, Δ + ℎ +1∕2, cos +1∕2 +1, - +1∕2 , Δ = ′ +1 , +1∕2 -, +1∕2 Δ + ℎ , +1∕2 ⎛ ⎜ ⎜ ⎝ +1∕2 , +1 - +1∕2 , Δ ⎞ ⎟ ⎟ ⎠ = - ′ ′ = 1 4 +1, +1∕2 + +1, -1∕2 + , +1∕2 + , -1∕2 ′ = 1 4 +1∕2, +1 + +1∕2, + -1∕2, +1 + -1∕2, (2)
where Δ , Δ define the the grid size on the longitude and latitude directions respectively; Δ is the time step; and

, , and indicate the values of their respective variable at time Δ and location ( Δ , Δ ). As shown by [START_REF] Brodtkorb | Efficient shallow water simulations on gpus: Implementation, visualization, verification, and validation[END_REF], this numerical scheme is well-suited for running in the GPU.

Given a domain and its bathymetry ℎ, the input required to integrate equations ( 1) is contained in the initial and boundary conditions. As mentioned by [START_REF] Wang | User manual for comcot version 1.7[END_REF], numerical stability is ensured by selecting Δ = × Δ ∕ , with < 1, where is the Courant-Friedrichs-Lewy number, = max √ ℎ and Δ is the minimum geodesic distance between two vertically or horizontally adjacent grid nodes.

Initial conditions

Although arbitrary initial conditions can be provided through ℎ, , and , there are formulas that characterize the generation of the tsunami depending on a set of parameters. Two of these are implemented in BROWNI: earthquakes and asteroids. For earthquake generation, we use a finite fault model, wherein the water is assumed to be at rest and a perturbation is added around the epicenter according to the formulas of [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. These formulas are an explicit solution to a linear elasticity problem that assumes that the area affected by the earthquake is rectangular and that it depends on parameters such as its length, width, hypocenter, fault slip, and the 3D orientation of the fault plane. Complex earthquakes can then be represented by superposing several of these rectangles with different parameters depending on the heterogeneity of the earthquake being modeled.

For asteroid-generated tsunamis, [START_REF] Ward | Asteroid impact tsunami: a probabilistic hazard assessment[END_REF] described several formulas that approximate the deformation of the water due to the impact of the asteroid depending on how its energy is converted into water waves. In BROWNI, one such formula is implemented. It also assumes that the initial surface has radial symmetry and can be described by a paraboloid of positive curvature whose size depends on the size, mass, density, and speed of the impactor, as well as the amount of energy that is transmitted to the water from the asteroid.

Boundary conditions

Three types of boundaries are considered: periodic, open, and closed. Periodic boundaries repeat the information of the opposite boundary and are used for the E-W borders whenever the domain covers a full circle of latitude. Open boundaries allow waves to leave the domain without reflections and are used whenever the domain does not cover a latitudinal circle. This boundary condition considers the trajectories of the characteristics of the Riemann invariants of the wave equations to extrapolate values in time, as explained by [27]. Closed boundaries are used to simulate a reflecting wall by imposing = = = 0; this is used to model a continental shoreline as a closed internal boundary, which is a reasonable assumption since, at a large scale, run-up displacements are negligible [27].

Software implementation

Data sources

Update the state of

Overview

BROWNI is designed to be integrated in web applications directly in the browser, as shown in Figure 1. First a user, such as an emergency manager or any person interested in visualizing and configuring tsunami simulations, interacts with a web application through a graphical user interface, on which relevant information is displayed, possibly along with interaction controls such as buttons, menus, and so on. Then, this application uses BROWNI to simulate tsunamis with the user configuration and application data, which is used to run the GPU instructions in the WebGL engine that runs the simulation efficiently. At the same time, state updates may come from other data sources, possibly producing changes in the simulator. What is important is that BROWNI is agnostic of the web-application's purpose or structure, as also of when user-driven or data-driven changes to the simulation are applied; these occur asynchronously, i.e., independently of the main flow of the simulation, and without requiring the user to leave the web browser.

BROWNI components details

To explain how these features are achieved in BROWNI, a "zoomed" architecture diagram is shown in Figure 2, detailing the components of BROWNI that facilitate the configuration, interactivity, and integration of the simulations.

First, the application should provide input data such as domain configuration, bathymetry, and initial conditions. Some of these data may be in different formats, such as CSV, JSON text files, or PNG and JPG images. This data is internally received by a reader component, which contains methods to parse these different formats. The reader converts this data into useful input for the shallow water model component, which can only receive input data in one format. Once the data is received, the shallow water model component is then in charge of executing the simulation sequentially. It then calls the WebGL instructions that use GPU acceleration to speed up the calculations.

To run several iterations until a specific timestamp is reached, a simulation controller component is implemented, which not only runs the simulation program, but also calls model functions to drive the simulation and extract results

to fill an HTML5 canvas element with a pseudo-color map or colored texture, representing the values of a desired variable such as current wave heights, maximum heights, or arrival times. Then, using the controller, the application can interrupt the simulation with common playback commands (play, pause, restart, etc. ) on user demand; and with the model, the application can collect the results of the simulation and populate its views, for example, by overlaying the canvas texture on a map or displaying a time series plot of a specific point of interest.

Asynchronous event-handling with life-cycle callbacks

Although the shallow water model and the simulation controller components help with interactivity and integration, it is important to note that they are not sufficient for handling asynchronous state changes in a predictable way. For example, one may be calling model functions before the data is available, or before the components exist, or it may be necessary to interrupt the simulation cycle under a certain specific condition, say, once the wave height is greater than a threshold or after a certain number of iterations. For this reason, we introduce the concept of "life cycle" shown in figure 3, that describes the sequence of relevant stages since the data is received until the controller has reached its final simulation time . The application can then subscribe to events by providing a set of callback functions (listed on the left side of figure 3) that are internally called at specific steps of the life cycle.

Model implementation in WebGL

The shallow water model component is the software representation of the algorithm and mathematical equations presented in section 2. It uses WebGL, a web standard for GPU accelerated computer graphics on the web browser [START_REF] Marrin | Webgl specification[END_REF], and as such, it requires one to adapt the code implementation to concepts and design constraints that are very specific to computer graphics. The chosen model is well suited for this approach as every node in the numerical mesh depends only on a small number of neighboring nodes, making it analogous to image rasterization algorithms [START_REF] Brodtkorb | Efficient shallow water simulations on gpus: Implementation, visualization, verification, and validation[END_REF].

The simulation (or "rendering") process is performed in a special type of program called Fragment Shader, which is written in a strongly-typed language called OpenGL Shading Language (GLSL), whose most basic purpose is to fill an RGBA matrix of pixels with the desired colors to build the target image, that is, a matrix with three intensity levels for red, green, and blue and one for the transparency level (alpha). In addition, by design, a Shader program has no memory of previously rendered scenes, therefore any prior information must be provided explicitly. For this reason, and for better performance, instead of rendering to the screen, BROWNI renders to a frame buffer object (FBO), which is a WebGL object that allows storing of the pixel data of a Shader into a texture that can be used as input later in other Shader instances. In total, it uses two FBOs:

1 to store the results of the previous time step and 2 for storing new results that depend on 1 . On each pixel of an FBO, the RGBA values are stored corresponding to ( , , , ℎ) at each point in the simulation.

Validation

The mathematical model presented in section 2 has been already validated by other authors, for example in [START_REF] Titov | Implementation and testing of the method of splitting tsunami (most) model[END_REF][START_REF] Liu | Numerical simulations of the 1960 chilean tsunami propagation and inundation at hilo, hawaii[END_REF] and more recently in [START_REF] Adams | Geoclaw model tsunamis compared to tide gauge results final report[END_REF][START_REF] Greenslade | An assessment of the diversity in scenario-based tsunami forecasts for the indian ocean[END_REF]. Here we examine the accuracy of the results against measurements of two different farfield tsunamis to ensure the correctness of the implemented model. This is also useful to show how BROWNI can be configured in a manual workflow where post-processing tasks are made outside of the browser context where the simulation runs.

To configure the simulations, it is necessary to characterize the scenario, provide output options, and optional life cycle callbacks to catch key events of the simulation. An example of a standalone HTML/JavaScript code that can be run from the web browser is shown in Figure 5. Two scenarios were studied, corresponding to the 8.8 Mw 2010 Maule earthquake and 9.1 Mw 2011 Tohoku earthquake. For each scenario, bathymetry and earthquake information were provided through external uncompressed binary files. The bathymetry source is ETOPO-1 [START_REF] Amante | Etopo1 1 arc-minute global relief model: procedures, data sources and analysis[END_REF], and the earthquake is represented by finite fault models proposed by [START_REF] Delouis | Slip distribution of the february 27, 2010 mw= 8.8 maule earthquake, central chile, from static and high-rate gps, insar, and broadband teleseismic data[END_REF] and [START_REF]M 9.1 -near the east coast of Honshu, Japan[END_REF] for the scenarios of Chile and Japan, respectively, using the formulas mentioned in section 2. As shown in Figure 5 Since the postprocessing is made outside the context where the simulation is running, the only provided lifecycle callback is the simulationDidFinish function (see figure 3), where it instructs the Controller to export files containing gridded values of maximum amplitude and arrival times after the simulation has reached the stopTime.

"Heat maps" of maximum amplitudes and travel time isochrones were produced in Python using the results of BROWNI. Figures 6 and8 show these results, which demonstrate how the implemented numerical model represents tsunami propagation, including reflection and refraction patterns that are developed due to its interaction with bathymetry and shore lines.

Line plots showing the time series of wave amplitude for each scenario are shown in Figures 7 and9. Black lines correspond to the measurements of Deep-ocean Assessment and Reporting of Tsunamis (DART) buoys after filtering out the astronomical tide using a low pass filter. DART buoys locations and identification numbers are also shown with circles in Figures 6 and8. The other lines represent results computed by different tsunami-simulation software: Easywave (red), a software developed in the German Research Center for Geosciences, Potsdam [START_REF] Christgau | A comparison of cuda and openacc: accelerating the tsunami simulation easywave[END_REF] that implements the same numerical scheme described in section 2 with C++ and CUDA; and GeoClaw (green), part of the ClawPack software [START_REF]Clawpack software[END_REF], that solves the non-linear shallow water equations using a finite-volume scheme designed to represent the propagation of tsunamis both at the global and regional scales, including inundation. All models' results were shifted in time to match the location of the first observed peak. The computed time-shift of each model with respect to the measurements is shown in tables 1 and 2. A negative time-shift indicates an early arrival.

In both scenarios, BROWNI and EasyWave show almost no difference from each other, since both solve the same set of equations using the same numerical scheme. GeoClaw shows an earlier arrival than BROWNI and Easywave, a behavior already observed in [START_REF] Adams | Geoclaw model tsunamis compared to tide gauge results final report[END_REF], that is explained by the different dispersion characteristics of the numerical schemes and equations. Finally, a consistent delay on the arrival of the first wave, that increases with the distance to the earthquake source, can be observed in all computed results. As reported by [START_REF] Abdolali | Role of compressibility on tsunami propagation[END_REF] and [START_REF] Watada | Traveltime delay and initial phase reversal of distant tsunamis coupled with the self-gravitating elastic earth[END_REF], this delay can be reduced by including additional physical features ignored in the fundamental assumptions of shallow water models such as water compressibility and earth's elasticity.

These results show that by using WebGL from the web browser as a computing engine, BROWNI can produce results as accurate as other tsunami simulation software packages. The next section exemplifies how additional versatility can be obtained, by taking advantage of WebGL as a web standard to produce interactive visualizations that are integrated with other sources of data and user-interface elements.

Application to scientific outreach and communication

We used BROWNI to try different kinds of interactive applications on outreach activities that aim at increasing tsunami risk awareness. For example, we have participated in the international competition "Mathematics of Planet Earth 2017" organized by Imaginary [START_REF]Winners of the second mathematics of planet earth competition[END_REF] and the Futur.E.S. festivals organized by Cap Digital [START_REF] Devillard | Les innovations à ne pas rater au festival futur.e.s en seine, ???[END_REF] (see picture (a) in figure 13). Specific details on outreach activities that we have performed are reported in [START_REF] Zamora | The 1730 great metropolitan chile earthquake and tsunami commemoration: Joint efforts to increase the country's awareness[END_REF].

By wrapping the simulator in a library we have separated concerns in what respects to tsunami-simulations and application development, thus simplifying the development process, making it easier and more sustainable. And once those applications were implemented, the most remarkable benefit has been the ability to reduce the total cost of gathering evidence from tsunami simulations to verify or discard previous assumptions. This makes it possible to focus the conversations on the scientific questions that are important for both the science communicators and the spectators instead of the complexities and implementation details of the simulator and the interactive visualizations.

The result is a dynamic and engaging experience for all participants. In this section we describe two of these applications that were developed with BROWNI, namely TsunamiLab, a Single-Page Application deployed at www.tsunamilab.cl, and the TsunamiLab-Pool, an interactive augmented reality device meant for public exhibitions such as science fairs.

TsunamiLab: a Single-Page Application

Tsunamilab, available at www.tsunamilab.cl, is a Single-Page Application (SPA), i.e., a web application loaded from a single HTML document and whose content is updated using Javascript. TsunamiLab possesses several features that enable users to observe the propagation of tsunami waves around the world, and select and modify scenarios based on historical and "synthetic" data, i.e., user-defined earthquake locations and magnitudes.

Figure 10 shows the main view of TsunamiLab, available at www.tsunamilab.cl. Besides BROWNI, Tsunami-Lab is built using three main tools: ReactJS [START_REF] Fedosejev | React. js Essentials[END_REF], a Javascript library maintained by Facebook that facilitates the development of complex interfaces under a reactive programming paradigm; Three.js [START_REF] Cabello | Three. js[END_REF] a popular general purpose 3D library that uses WebGL; and the USGS Earthquake Catalog web API, an API that allows custom searches for earthquake information using a variety of parameters [START_REF] Survey | API Documentation -Earthquake Catalog[END_REF].

Frontend architecture with BROWNI

The software architecture is described in a component tree as shown in figure 11. On this tree, components with a common parent can access its state variables and pass them to the next component in the same branch. Components down the tree on a same branch can also receive callback functions to trigger changes in state variables of a parent component, triggering updates and data requests. For example, the magnitude and location of the current scenario are shared state-variables between the BROWNI, canvas, scene, synth. scenarios, and historical scenarios components, and thus any time there is an update (triggered by an user interaction, for example) all of these components will be notified and re-rendered.

TsunamiLab uses BROWNI to run the simulation and Three.js for displaying a 3D scene (in the scene component of Fig. 11) where the globe, pins, and other 3D graphics components are rendered. The integration of BROWNI with Three.js is performed by sharing an HTMLCanvasElement (the canvas component of Fig. 11) between them.

This canvas serves for creating instances of the WebGL context used by BROWNI to run the simulations, and also for storing the "heatmap" with the colors of the simulation at every timestep, which is then displayed in the globe of the scene as a texture on top of the base map.

Numerical model settings

The

Historical scenarios

Historical scenarios are represented in the globe as a collection of wave pins scattered around the globe at earthquake locations. Earthquake location, magnitude, and focal mechanism information are obtained through the publicly available USGS Earthquake Catalog web API. Whenever the user clicks on a visible pin, an HTTP request is sent to the USGS Earthquake Catalog and the shared magnitude and location states are updated, triggering the update of BROWNI with the new scenario. At the same time, text information is updated in the "Scenario Info" component (bottom right table in figure 10) to give the user context about the current scenario.

Synthetic scenarios

Synthetic scenarios can be configured in a separate box in the top right corner of figure 10. The magnitude can be chosen by increasing or decreasing the radius of the inner circle, in which case the magnitude of the new scenario is selected to be proportional to 3 . Earthquake location can also be selected by dragging and dropping the orange pin on top of the globe. BROWNI is updated automatically whenever the location or the magnitude are changed through either of these two actions.

TsunamiLab-Pool: An augmented reality-device

With the development of TsunamiLab we noticed that people showed a greater interest after trying out their own scenarios to verify or discard their previous assumptions about tsunamis. However, from our experience in various exhibitions, we found that the classic mouse-keyboard-monitor system for visualization and interaction was still uncomfortable to use in contexts such as science fairs or exhibitions at local schools.

For this reason we created the TsunamiLab-Pool, an augmented reality device in which the visualization of the simulation is projected onto a circular surface supported by a rigid frame. The interaction is handled by an external controller, to rotate the globe and change the location and magnitude of new earthquakes. We tried two different controllers: the Sony PlayStation 3 Move (PSMove) controller and the Leap Motion controller. The former is a wireless controller developed by Sony that, in addition to traditional buttons, includes a marker whose position is tracked by a camera and mapped accordingly into the scene to help with pointing actions. The latter consists of a set of infrared lights and cameras that are used to infer the position of the user's fingers, hands, and arms. With the Leap Motion Controller the user is actually able to control the parameters of the simulation and the visualization by just moving the hands in the air (see picture (b) of figure 13).

The TsunamiLab-Pool software is built as in diagram 12. First the user interacts with the physical device where the simulation is projected using a controller. This controller sends the raw information of the interaction to a local controller-service. This controller-service then translates the data from the raw interactions into simulation commands that are forwarded to the simulation web app, which in turn, uses them to configure and calculate the new scenario and render the visualization back into the projection of the physical device.

The controller service is a python program that extracts the data from the actual remote controllers. It uses the PSMoveAPI python library [START_REF] Perl | Cross-platform tracking of a 6dof motion controller, Using Computer Vision and Sensor Fusion[END_REF] for the PSMove controller and the python API of the Leap Motion SDK version 2.3, for the Leap Motion controller. The simulation web app is a plain HTML-Javascript-CSS application, served locally with a simple http server that receives messages from the controller-service in real-time through a websockets connection.

These messages change the state of the simulation, that is running with BROWNI, and also the visualization, which is rendered using Three.js.

Thanks to the API of BROWNI, no changes to the source code of the library were necessary. This made the development experience cleaner and more sustainable than, for example, duplicating the source code on each of the applications, since it allows for a better separation of concerns. And because BROWNI uses WebGL, the integration inside of the website was straight-forward and its on-site deployment only required a gaming computer to give a smooth experience to visitors and science communicators.

Discussion on Software and Hardware Limitations

Running the tsunami model on the browser using WebGL offers several opportunities as disussed before, but some limitations or constraints should be discussed.

Multi-threading

Using multiple threads or processes can be beneficial for a webpage to prevent different tasks from interrupting each other inside the browser. This is possible in recent versions of web-browsers such as those based on Chromium, by using two JavasScript API's: the WebWorker API to run tasks in the background, detached from the main thread; and the OffScreenCanvas, to run the WebGL simulation from inside a WebWorker, instead of an HTMLCanvasElement.

Naturally, there are two ways on which multi-threading can improve the execution of the simulation: preventing other tasks from slowing down the computations; and preventing the simulation from slowing down the other tasks.

Although the first case is straight-forward, and is the main use-case for the WebWorker API, the second case must be reviewed carefully when designing a web application, since large simulations may still interrupt the webpage and other applications outside of the browser. This happens because GPUs are designed to handle only one rendering task at a time, including the rendering of the display's content. This does not happen with CPUs since operating systems usually employ preemptive multitasking, a feature that allows to temporarily pause heavy tasks while keeping the computer responsive. This feature is not usually available in GPUs [START_REF] Tanasic | Enabling preemptive multiprogramming on gpus[END_REF] and the operating system may even restart the GPU driver, assuming that the GPU was blocked, even though it was only processing a heavy task [START_REF]Microsoft, Resource Limits (Direct3D 11) -Windows app developer documentation[END_REF]. The only workaround in WebGL is to make a compromise between the size of the simulation (mesh size or number of time-steps computed per second) and the responsiveness of the system.

To overcome this limitation, native APIs may make use of a secondary GPU if available, however, this is not possible in WebGL since the browser can only use the same GPU as the user's display.

Background tabs and batch simulations

Leaving a simulation running in a background tab may cause the browser to temporarily slow it to even less than one frame per second, in order to alleviate resources for the currently active browsing tab. For cases where interactive simulation and visualization is the priority, as is the interest on this work, this is convenient since it prevents the user from missing important changes while they finish other tasks in the new tab. However, this can be limiting in other cases, for example, if one seeks to run batch simulations in the background, and process their results only at the end.

In this case, users could disable this feature using the command-line options of the browser's executable, which is the case for Chrome since version 57 [START_REF] Tanasic | Enabling preemptive multiprogramming on gpus[END_REF], for example. Disabling this feature will thus compromise the setup-cost and browser support, so its convenience should be evaluated on a case-by-case basis.

A different approach for running batch simulations is to, instead of running the simulations in the background, run them interactively in the forefront while changing their parameters on the fly. For smaller scenarios, several instances of BROWNI could be created in the same browser, and for more complex ones, a tiled-display visualization could be built similarly to [START_REF] Hsieh | Parallel tsunami simulation and visualization on tiled display wall using opengl shading language[END_REF] and [START_REF] Kamakshidasan | Comparative visualization of deep water asteroid impacts on ultrahigh-resolution wall displays with seawall[END_REF].

Scalability

Using CUDA to scale-up the size of the simulation has already been shown to be possible, for example in [START_REF] Macías | Comparison and computational performance of tsunami-hysea and most models for lantex 2013 scenario: Impact assessment on puerto rico coasts[END_REF][START_REF] De La Asunción | Numerical simulation of tsunamis generated by landslides on multiple gpus[END_REF].

However, in the case of WebGL this is not straight-forward, if ever possible or convenient.

In a single node, WebGL has some limitations, which are already listed in [START_REF]GFXFundamentals, WebGL2 Cross Platform Issues[END_REF]. These limitations also exist in the native APIs of OpenGL [START_REF] Group | Opengl® 4.5 reference pages[END_REF], OpenGL ES [START_REF] Group | Opengl® es 3.2, ???[END_REF], and Direct3D [START_REF]Microsoft, Resource Limits (Direct3D 11) -Windows app developer documentation[END_REF]. Each GPU has a maximum texture size that it can store and it depends not only on the available memory of the GPU but also on the specific design of the device.

Also, even if the information were split into several textures by decomposing the domain, there exists a maximum number of textures that the GPU can store that also depends on the particular device. These factors should be taken into consideration in order to give proper support to the platforms of each web application. Some additional developments could be done to overcome these limitations. For example, one could try to decompose the simulation among different computers and WebGL contexts with a real-time communication protocol such as WebSockets to step through the simulation, alleviating the resources of each computer. Also, adding server-side support could be done with the Nodejs JavaScript runtime instead of the browser's. However, these approaches are outside the scope of this work and have not been studied here.

Browser performance

Browsers such as Firefox, Edge, Google Chrome and others based on Chromium run WebGL calls using the Almost Native Graphics Layer Engine (ANGLE) API [START_REF] Project | ANGLE -Almost Native Graphics Layer Engine[END_REF], translating them into one of the graphics APIs available for that platform; usually OpenGL on Linux and Mac, and Direct3D on Windows. This naturally means that any WebGL application will be at most as fast as its native counterpart.

To examine the performance of BROWNI on the web browser, the elapsed time for the 8.8 Mw Maule, Chile, 2010 case (introduced in detail later on section 4) is shown in Table 3. Simulations were performed on a gaming laptop with an Intel Core i7-7700HQ 2.8 GHz CPU and 8 GB of RAM with Intel HD Graphics 630 integrated graphics card and also with an NVIDIA Geforce GTX 1060 dedicated graphics card. Two grid sizes were used with a spacing of 3 and 15 minutes until 25 hours of simulation passed. Although the execution times were relatively similar for the coarse grid, the simulation ran 5.16 times faster with the dedicated GPU. This speedup is smaller than, for example, the 10.16 speedup reported in [START_REF] Christgau | A comparison of cuda and openacc: accelerating the tsunami simulation easywave[END_REF] for EasyWave running on CUDA, instead of on the CPU, which could be explained by the previously mentioned overhead.

Although it can also be influenced by other factors, such as differences in the architectures of CPUs and integrated graphics cards, this still illustrates that speedup differences are expected in practice. However, as was described in the previous section, for applications in science education and risk communication, the interactivity and efficient visualization can increase the performance of other tasks needed to use the simulator, making the user more proficient on understanding its results.

Conclusions

On this study, it has been shown how GPU computing from web browsers can allow users to gather and analyze data from simulations efficiently, without being concerned by implementation details and tasks that are not relevant for understanding the physical phenomena of tsunami propagation and its potential implications to tsunami risk management, such as setup, configuration, pre-processing and post-processing of the results.

To demonstrate this, a Javascript library called BROWNI was introduced, which has an API designed to facilitate the creation of interactive GPU-powered tsunami simulations and visualizations in the web browser. By using this API, it was possible to separate concerns when building web applications, facilitating code-reusability and maintainability across different software platforms. Some limitations of this approach were also examined, finding that large scenarios may experience reduced performance and reduced mesh-resolution scalability support, as also necessary trade offs on browser-support to benefit from CPU multi-threading and background processing.

Results were also validated with other simulation software and measurements. These are coherent with the fact that BROWNI implements algorithms already used in operational and research settings: differences found on time-series were recalled, however these had already been discussed in the literature. This confirms the correct implementation of the algorithms.

Finally, we showed two software platforms that were developed using BROWNI, whose aim is to increase tsunami risk awareness: TsunamiLab and TsunamiLab-Pool. By creating and using these applications on different opportunities as tsunami-risk communication tools it was possible to confirm that: (1) that the benefits of the interactive simulation and visualization on the web browser can overcome its limitations, since users can focus on the scientific questions instead of the implementation details, thus producing a more engaging experience; (2) that developers can benefit from a more sustainable developer experience by using the BROWNI API as a part of their applications without worrying about the source code or other details of the simulator.

Future work may include changes to the model component to increase the accuracy in the estimation of arrival times, as proposed by [START_REF] Watada | Traveltime delay and initial phase reversal of distant tsunamis coupled with the self-gravitating elastic earth[END_REF]. In addition, improving the numerical dispersion to match the physical dispersion of linear waves in the ocean as proposed by [START_REF] Ha | Tsunami propagation over varying water depths[END_REF] and ultimately the inclusion of higher resolution algorithms to include nonlinear effects and bottom friction is important to represent smaller scale phenomena. Other questions that remain open are in regards to the possibility of using WebGL for high performance computing in tsunami modeling, which would require adding server side support and distributed computing capabilities to extend the scalability of the simulator as discussed on section 6. Furthermore, the question of the effectiveness of interactive tsunami simulations to support educational and/or communicational methodologies for improving tsunami awareness, risk perception, and other topics has not been covered here and could also be studied with this library.

Figure 1 :

 1 Figure 1: BROWNI architecture overview for a standard use case

Figure 2 :

 2 Figure 2: Architecture of BROWNI in terms of its internal components

Figure 4

 4 Figure 4 shows the flow of the simulation, in the scope of the shallow water model component. After the data is received from the reader component (top diagram), the bathymetry is rendered into a texture, the simulation timestep is calculated and the initial condition is rendered into 1 . Then, when the runSimulationStep function of the shallow water model component is called, 1 is assigned as the previous step texture and the simulation shader program is run, with its output stored into 2 . The flow ends with 1 and 2 being swapped, such that the data is ready when the runSimulationStep function is called again.

Figure 3 :Figure 4 :

 34 Figure 3: Life cycle diagram of a simulation in BROWNI.

Figure 5 :

 5 Figure 5: Example standalone HTML file with Javascript code for configuring BROWNI to produce heat maps of figures 6 and 8

  , the domain covers the spherical rectangle [90, 325.83] × [-60, 70] and is discretized in a spherical uniform grid of 4717 × 2600 ≈ 12.26 million nodes with a spacing of 3 min. Each simulation is run until 25 h of propagation (the stopTime parameter) with a timestep of 2.9365 s, which is configured by default for a CFL number of 0.5 to respect the stability condition explained in section 2; all boundary conditions are open. Computed results are shown in Figures 6, 7, 8, and 9, which correspond to the 8.8 Mw 2010 Maule earthquake and 9.1 Mw 2011 Tohoku earthquake.

Figure 6 :

 6 Figure 6: "Heat maps" of maximum wave amplitude and arrival time isochrones calculated by BROWNI for the 8.8 Mw scenario of Maule, Chile, 2010 and DART buoys locations.

Figure 7 :

 7 Figure 7: Time series of wave amplitude for the 8.8 Mw scenario of Maule, Chile, 2010, obtained by BROWNI (blue line), EasyWave (red line), GeoClaw (green line) and de-tided measurements (black line) of DART buoys. Time-shifts applied are shown in table 1

Figure 8 :

 8 Figure 8: "Heat maps" of maximum wave amplitude and arrival time isochrones calculated by BROWNI for the 9.1 Mw scenario of Tohoku, Japan, 2011 and DART buoys locations.

Figure 9 :

 9 Figure 9: Time series of wave amplitude for the 9.1 Mw scenario of Tohoku, Japan, 2011, obtained by BROWNI (blue line), EasyWave (red line), GeoClaw (green line) and de-tided measurements (black line) of DART buoys. Time-shifts applied are shown in table 2

  simulation domain covers the spherical rectangle ( , ) ∈ [-180, 180] ×[-70, 70] with a numerical grid of 10 minutes resolution. Bathymetry is sampled from the ETOPO-1 dataset and compressed into a PNG image with the same resolution as the simulation to facilitate navigation. The quantization error from the PNG compression is found to be insignificant, since the only region where it becomes important is at the shores, which are excluded from the numerical domain of the simulator. The E-W boundary conditions are periodic (at = -180, 180), whereas the N-S borders are open (at = -70, 80); as described on section 2, continental shorelines are simulated as inner reflecting walls.

Figure 10 :Figure 11 :

 1011 Figure 10: Screenshot of the main view of www.tsunamilab.cl

Figure 12 :

 12 Figure 12: Diagram with the architecture of the TsunamiLab-Pool.

Figure 13 :

 13 Figure 13: Pictures of activities performed with non-experts users: (a) Interactive display that uses the SONY ™PS Move Controller to configure earthquake magnitude and location, at the 2018 Paris FUTUR.E.S. technology festival organized by Cap Digital; (b) Interactive display that uses the Leap Motion controller to configure earthquake magnitude and location using hand gestures at a local science fair in Chile.

Table 1

 1 Time-shift applied to the results of each simulation-software of figure7sorted by arrival-time at each buoy (column "Arrival"). A negative amount indicates an earlier arrival.

			44	0.48	-1.69
	43412	9.81	-2.06	-1.35	-5.85
	46412	13.15	-2.87	-2.46	-8.96
	51425	14.43	-13.20	-13.21	-18.14
	46411	14.48	-10.06	-9.59	-15.98
	46407	14.98	-4.05	-3.58	-11.04
	46404	15.60	-6.57	-6.42	-13.93
	21413	21.73	-23.83	-19.41	-27.24
	52403	22.11	-12.14	-11.17	-15.55

table 1

 1

		Arrival	BROWNI EasyWave	GeoClaw
	Buoy ID (hours) (minutes) (minutes) (minutes)
	21413	1.34	-1.90	-2.00	-2.57
	52403	5.18	-6.68	-6.22	0.53
	51425	8.59	-4.64	-4.57	-8.25
	46404	8.89	-8.33	-8.19	-10.35
	46407	8.99	-9.49	-9.10	-10.71
	46411	9.34	-9.83	-10.05	-13.48
	46412	10.36	-8.81	-8.69	-11.83
	43412	13.50	-6.91	-6.89	-11.97
	32412	19.45	-13.81	-12.83	-18.37

Table 2 Time

 2

-shift applied to the results of each simulation-software of figure

9

sorted by arrival-time at each buoy (column "Arrival"). A negative amount indicates an earlier arrival.

Table 3

 3 Execution times for the simulation of the 8.8 Mw, 2010 Maule earthquake for two grids of different resolution, using an integrated and dedicated GPU. Integrated GPU: INTEL(R) HD Graphics 630; Dedicated GPU: NVIDIA Geforce GTX 1060.

		Integrated GPU Dedicated GPU
	15 minutes 943 × 520 grid	1.6 min	1.44 min
	3 minutes 4717 × 2600 grid 43.6 min	8.44 min

Acknowledgements

This project has been funded by projects CONICYT/FONDAP/15110017 "Centro de Investigación para la Gestión Integrada de Desastres Naturales" and CORFO 10CEII-9157 "Inria Chile". Support was also received by the Marine Energy Research & Innovation Center (MERIC, project CORFO 14CEI2-28228).

The authors are also grateful for the contributions of professor Felipe Cortez and the students of the course "Usabilidad y Nuevos Medios" of the School of Design of the second semester of 2018 of the Pontifical Catholic University of Chile, for their feedback and design proposals for TsunamiLab.

Software availability of BROWNI

Software availability

The source code of BROWNI is available under an open-source MIT license at https://github.com/jgalazm/ browni, and detailed information is specified in Table 4. Visualizations and data processing were performed in Python with matplotlib, and the source code for the cases presented herein is written in Jupyter Notebooks and versioned at https://github.com/jgalazm/comp-geo-paper-figures.


Mike J Newland
email: m.newland@uea.ac.uk

Patricia Martinerie

Emmanuel Witrant

Detlev Helmig

David R Worton

Chris Hogan

William T Sturges

Claire E Reeves

Changes to the chemical state of the northern hemisphere atmosphere during the second half of the twentieth century

The NO X (NO and NO 2 ) and HO X (OH and HO 2 ) budgets of the atmosphere exert a major influence on atmospheric composition, controlling removal of primary pollutants and formation of a wide range of secondary products, including ozone, that can influence human health and climate. However, there remain large uncertainties in the changes to these budgets over recent decades. Due to their short atmospheric lifetimes, NO X and HO X are highly variable in space and time, and so the measurements of these species are of very limited value for examining long term, large scale changes to their budgets. Here, we take an alternative approach by examining long-term atmospheric trends of alkyl nitrates, the formation of which is dependent on the atmospheric NO/HO 2 ratio. We derive long term trends in the alkyl nitrates from measurements in firn air from the NEEM site, Greenland.

Their mixing ratios increased by a factor of 4 -5 between the 1970s and 1990s. This was followed by a steep decline to the sampling date of 2008. Moreover, we examine how the trends in the alkyl nitrates compare to similarly derived trends in their parent alkanes (i.e. the alkanes which, when oxidised in the presence of NO X, lead to the formation of the alkyl nitrates). The ratios of the alkyl nitrates to their parent alkanes increase from

Introduction

The NO X (NO + NO 2 ) and HO X (OH + HO 2 ) budgets of the troposphere act to control the concentrations of oxidants such as OH, ozone and NO 3 (Fig. 1). These in turn control removal of pollutants from the atmosphere. Emissions of NO x in the northern hemisphere are mainly anthropogenic, with roughly equal proportions from power generation and transport [START_REF] Olivier | Global emissions sources and sinks[END_REF][START_REF] Olivier | Global emissions sources and sinks[END_REF]. NO X and HO X are linked through the production of ozone, which is positively correlated with NO X concentrations in the background atmosphere through the photolysis of NO 2 (Reactions R1-R2). The photolysis of ozone in the presence of water vapour then leads to the production of OH (Reactions R3-R4). Removal of NO X from the atmosphere is controlled by the reaction of NO 2 with OH during the daytime (Reaction R5). This forms nitric acid, HNO 3 , which is lost from the atmosphere by wet deposition. At night, and during the winter, the heterogeneous reaction of the NOx reservoir species N 2 O 5 (formed from the reaction of NO 2 with NO 3 (Reactions R6-R7) with H 2 O on aerosol becomes an important NO X sink (Reaction R8). OH and HO 2 rapidly interconvert through the reactions of OH with CO and hydrocarbons, such as alkanes, and the reaction of HO 2 with NO (Reaction R9). The reaction of NO with peroxy radicals (HO 2 and RO 2 -Reactions R9-R10) recycles the NO X back to NO 2 . The main removal processes for HOx are the reaction of OH with NO 2 (Reaction R9) and the reaction of HO 2 with peroxy radicals (including self-reaction) (Reactions R11-R12). However, changes to the atmospheric concentrations of both HO X and NO X during the previous century are poorly constrained. This is because all HO X and NO X species are short lived, present at low concentrations (0.01 -10 ppt) and have a high spatial and temporal variability. (e.g. [START_REF] Stone | Tropospheric OH and HO 2 radicals: field measurements and model comparisons[END_REF]. This makes them difficult to measure and trends difficult to identify (based on spatially and temporally variable data sets). Furthermore, a range of state of the art atmospheric chemistry transport models give no consensus of even the sign of OH change between 1850 and 2000 [START_REF] Naik | Preindustrial to present-day changes in tropospheric hydroxyl radical and methane lifetime from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP)[END_REF]. However, the models do agree that between 1980

HO X sources

NO 2 à NO + O (R1) O + O 2 à O 3 (R2) O 3 à O( 1 D) + O 2 (R3) O( 1 D) + H 2 O à 2OH (R4)
and 2000 there has been an increase in northern hemisphere OH concentrations, with the best estimate of the increase being 4.6 (± 1.9) %. This modelled increase is driven by increases in the NO x burden and in the water vapour concentration.

To attempt to study historical trends in HO X and NO X we have examined trends in longer living species which are affected by changes to HO X and NO X in the atmosphere.

In this paper we report long term atmospheric trends of three alkyl nitrates derived from Arctic firn air. These are chemically produced in the atmosphere from the oxidation of alkanes and subsequent reaction of the peroxy radical formed with NO. The alkyl nitrate records are combined with previously reported trends of their parent alkanes from the same Arctic firn site.

These records provide a proxy from which we can learn about the chemical state of the atmosphere in which they were formed.

Alkanes

Emissions of butanes (C 4 H 10 ) and pentanes (C 5 H 12 ) to the atmosphere are almost entirely anthropogenic (> 98 % globally [START_REF] Pozzer | Observed and simulated global distribution and budget of atmospheric C2-C5 alkanes[END_REF]), associated with fugitive emissions during oil and natural gas drilling and evaporation and combustion of fossil fuels such as in road vehicles [START_REF] Pozzer | Observed and simulated global distribution and budget of atmospheric C2-C5 alkanes[END_REF][START_REF] Pétron | Hydrocarbon emissions characterization in the Colorado Front Range: A pilot study[END_REF]Helmig et al., 2014a). The primary removal mechanism of these species from the atmosphere is reaction with the hydroxyl radical, OH (minor sinks include reaction with atomic chlorine, Cl, and the nitrate radical, NO 3 ).

Atmospheric mixing ratios of butanes and pentanes display a large seasonal cycle in mid-high latitudes (e.g. [START_REF] Swanson | Seasonal variation of C2-C4 nonmethane hydrocarbons and C1-C4 alkyl nitrates at the Summit research station in Greenland[END_REF][START_REF] Helmig | Volatile Organic Compounds in the Global Atmosphere[END_REF] due to changes in their chemical lifetimes (~1 month in the winter and 4 -5 days in the summer) driven by the seasonal cycle in OH concentration. Emissions are not thought to have a significant seasonal variability [START_REF] Pozzer | Observed and simulated global distribution and budget of atmospheric C2-C5 alkanes[END_REF].

Alkyl Nitrates

Alkyl nitrates (RONO 2 ) are secondary oxidation products of the alkanes (RH). Their atmospheric lifetimes are on the order of months in winter and ten days in summer (Clemitshaw Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3. 0 License. et al., 1997). Consequently, they display a strong seasonal cycle in the Arctic, with peaks in the late winter/early spring and minima in the summer [START_REF] Swanson | Seasonal variation of C2-C4 nonmethane hydrocarbons and C1-C4 alkyl nitrates at the Summit research station in Greenland[END_REF], similar to the alkanes. Alkyl nitrates are formed when alkanes react with OH to form a peroxy radical, RO 2 (Reaction R13), which subsequently reacts with NO to form an alkyl nitrate (Reaction 14b) (e.g. [START_REF] Talukdar | Atmospheric fate of several alkyl nitrates. 1. Rate coefficients of the reactions alkyl nitrates with isotopically labelled hydroxyl radicals[END_REF]. The '2008 EU hole' was sampled using the firn air system of the University of Bern [START_REF] Schwander | The age of the air in the firn and the ice at Summit, Greenland[END_REF], and the 'US' hole, sampled using the US firn air system [START_REF] Battle | Atmospheric gas concentrations over the past century measured in air from firn at the South Pole[END_REF]

Firn Analysis

The firn air samples from the 'EU' hole at NEEM were analysed for alkyl nitrates at UEA using a GC-MS in Negative Ion Chemical Ionisation mode (GC-NICI-MS) (e.g. [START_REF] Worton | Gas chromatography negative ion chemical ionization mass spectrometry: application to the detection of alkyl nitrates and halocarbons in the atmosphere[END_REF].

2-pentyl nitrate and 3-pentyl nitrate are presented together as 2+3-pentyl nitrate because the two are not baseline separated in the chromatogram.

The NEEM samples were analysed using the UEA calibration scale. This was converted to the NCAR scale (against which the North GRIP 2-butyl nitrate and 2+3-pentyl nitrate samples are calibrated) for direct comparison with the North GRIP atmospheric histories from Worton et al.

(2012) and with the in-situ measurements at Summit by UCI [START_REF] Swanson | Seasonal variation of C2-C4 nonmethane hydrocarbons and C1-C4 alkyl nitrates at the Summit research station in Greenland[END_REF][START_REF] Dibb | An overview of air-snow exchange at Summit, Greenland: Recent experiments and findings[END_REF]. This scaling was based on an inter-comparison between the UEA and NCAR standards in 2005 and 2012/13. These led to a rescaling of the UEA 2-butyl nitrate values by 1.245 and 2+3-pentyl nitrate by 1.409. The measurements of 3-methyl-2-butyl nitrate were not rescaled as the North GRIP measurements were made on the UEA scale.

Firn air samples from the 'EU' hole at NEEM were analysed for alkanes at the Max Planck Institute Laboratory (MPI) by gas chromatography with flame ionisation detection (GC-FID) (see [START_REF] Baker | Analysis of non-methane hydrocarbons in air samples collected aboard the CARIBIC passenger aircraft[END_REF] for further details). At the Institute of Arctic and Alpine Research (INSTAAR) firn air samples were analysed samples from both the 'EU' and 'US' holes at NEEM for alkanes by GC-FID (see [START_REF] Pollmann | Sampling, storage, and analysis of C2-C7 non-methane hydrocarbons from the US National Oceanic and Atmospheric Administration Cooperative Air Sampling Network glass flasks[END_REF]Helmig et al., 2014b for further details).

Firn Modeling

The air sampled from any given depth in the firn column is representative of a range of ages because of the inter-connected nature of the firn. Firn models can be used to derive the atmospheric history of a gas from measurements of air trapped in the firn. The extent and rate, at which the gas diffuses through the firn, depends on the diffusivity profile of the firn, the diffusivity coefficient of the gas, and on the gravitational fractionation (caused by the molecular weight) of the gas. The diffusion profile is different for every firn site. For determining the atmospheric history of a gas from firn air measurements, the firn diffusion profiles must first be constrained. This is done using a series of reference gases with well known atmospheric histories. At NEEM the reference gases used were CO 2 , CH 4 , SF 6 , HFC-134a, CFC-11, CFC-12, CFC-113, and CH 3 CCl 3 , as well as 14 CO 2 [START_REF] Witrant | A new multi-gas constrained model of trace gas non-homogeneous transport in firn: evaluation and behaviour at eleven polar sites[END_REF].

Each gas also has a different diffusion rate through the firn based on its molecular structure, this is called the diffusion coefficient. The diffusion coefficient is calculated relative to a reference gas, generally CO 2 . Different methods have been reported for the calculation of these diffusion coefficients (e.g. [START_REF] Chen | New generalized equation for gas diffusion coefficient[END_REF][START_REF] Fuller | A new method for prediction of binary gasphase diffusion coefficients[END_REF]. The diffusion coefficients of the alkyl nitrates were calculated using the method of [START_REF] Fuller | A new method for prediction of binary gasphase diffusion coefficients[END_REF] based on the sum of the Le Bas molar volumes of the molecule. The diffusion coefficients used for the firn modelling for each molecule within this work are given in Table 1.

The inverse model used for the atmospheric history reconstructions was the most recent version of the LGGE-GIPSA atmospheric trend reconstruction model described in [START_REF] Witrant | Input Estimation from Sparse Measurements in LPV Systems and Isotopic Ratios in Polar Firns[END_REF].

The atmospheric mole fraction derived from the firn reconstructions represents an annual mean.

The alkanes and alkyl nitrates examined in this work have a strong atmospheric seasonality due to changes in their chemical lifetimes driven by seasonal variability in OH concentration in the air masses in which they are transported to the Arctic. Thus the seasonal cycle of both species follows a roughly sinusoidal curve with a peak in the late winter (March) and a trough in midsummer (July-August) [START_REF] Swanson | Seasonal variation of C2-C4 nonmethane hydrocarbons and C1-C4 alkyl nitrates at the Summit research station in Greenland[END_REF]. Consequently, changes to the firn derived mole fractions are likely to be dominated by changes to winter-time atmospheric concentrations.

This model cannot take into account the seasonality in the signal that is preserved in the upper part of a firn profile. Therefore, measurements above a certain depth must be excluded from the Table 1 Diffusion coefficients used in the firn modeling, calculated from Le Bas molecular volumes using the method of [START_REF] Fuller | A new method for prediction of binary gasphase diffusion coefficients[END_REF].

Alkyl Nitrate Trends

Atmospheric histories of the three alkyl nitrates 2-butyl nitrate, 2+3-pentyl nitrate, and 3methyl-2-butyl nitrate (formed from n-butane, n-pentane, and iso-pentane, respectively) derived from firn air measurements at NEEM are shown in Figure 2. The records of all three alkyl nitrates show similar features (as would be expected from the similar sources and sinks).

All show a steep increase in mixing ratio from the 1970s to the 1990s with increases of a factor of 4 -5. The peak in the 1990s is followed by a steep decline to the sampling date of 2008. Considering Reactions R13 -R16, the trends in the alkyl nitrate mixing ratios (Fig. 2) could be caused by:

(i) Changes to the atmospheric mixing ratios of the parent alkanes;

(ii) Changes to [OH].t, where t is time since emission of the alkane. i.e. the amount of photochemical processing that the air mass in which the alkyl nitrates are being formed undergoes before reaching the Arctic;

(iii) Changes to the production efficiency of the alkyl nitrates, i.e. whether the peroxy radical reacts with NO (Reaction R14) or with HO 2 (Reaction R11;

(iv) Changes to the alkyl nitrate sinks, i.e., changes in [OH] or radiation.

Concerning point (i), the peak in alkyl nitrate mixing ratios in the 1990s is not contemporaneous with that of the parent alkanes (~1980 -Figure 3). This suggests that the changes to the alkyl nitrate mixing ratios are not being primarily driven by changes to the parent alkane. By considering the ratio of the alkyl nitrate to its parent hydrocarbon, using the firn derived alkane trends from NEEM presented in Helmig et al. (2014b), we can effectively remove the effect of changes to the parent hydrocarbon from the alkyl nitrate signal. This is done in Sections 4, 5 and 6, when points (ii) and (iii) are explored further. Concerning point (iv), there is evidence for global dimming (i.e. a decrease in surface solar radiation) of about 5 % between 1960 and 1990 in the northern hemisphere. However, this began to turn around during the mid 1980s and there was a brightening trend between 1985 and 2000 [START_REF] Wild | From dimming to brightening: decadal changes in solar radiation at earth's surface[END_REF]. This minor change to the alkyl nitrate sink is unlikely to have had a noticeable effect on mixing ratios.

Points (ii) and (iii) are discussed in the following sections. The trend shows similar features to those of the alkyl nitrate trends. The ratio increases from the 1970s to the late 1990s and then subsequently declines.

To investigate the drivers that might have led to these changes in [NO]/[HO 2 ] ratio, we shall now examine how the NO x and HO 2 concentrations may have changed.

Changes to atmospheric NO x concentrations

NO x sources

The atmospheric NO X concentration is determined by the relative magnitudes of the sources and sinks. The main sources of NO X in the northern hemisphere are anthropogenic emissions from fossil fuel use, power stations and transport [START_REF] Olivier | Global emissions sources and sinks[END_REF][START_REF] Olivier | Global emissions sources and sinks[END_REF]. Assuming that these bottom up emissions estimates are correct in the timing of the NO X emissions changes, it seems unlikely that the increase in the [NO]/[HO 2 ] ratio derived for the period 1970 to the late 1990s was driven primarily by changing NO X emissions.

The subsequent decline in [NO]/[HO 2 ], on the other hand, may well have been driven by decreasing emissions. Measurements at a range of UK sites showed a decrease in NO X concentrations from 1996 (the beginning of the reported measurements) to 2002 -2004, of 1 -3.5 % per year, depending on the site [START_REF] License | Trends in NO x and NO 2 emissions and ambient measurements in the UK[END_REF]. During the same period, the fraction of the NO X that is NO 2 (f-NO 2 ) roughly doubled. Recent trends at many European sites show similar trends with small decreases in NO x between 1999 and 2007 (the period for which measurements are available) but level or increasing NO 2 through the same period [START_REF] License | Trends in NO x and NO 2 emissions and ambient measurements in the UK[END_REF][START_REF] Gilge | Ozone, carbon monoxide and nitrogen oxides time series at four alpine GAW mountain stations in central Europe[END_REF] -suggesting that NO has decreased by more than NO x concentrations.

Declining NO X emissions have been used to explain these trends in measured concentrations.

However, the decline in these ambient concentrations is not as large as would be expected using current emission inventories [START_REF] License | Trends in NO x and NO 2 emissions and ambient measurements in the UK[END_REF].

NO x sinks

At mid-high latitudes, in the daytime, during the summer, the main sink for NO X is the reaction of NO 2 with OH. This reaction produces nitric acid (HNO 3 ), much of which is then removed from the atmosphere by wet deposition. However, at night and during the winter months, when daily mean [OH] is more than an order of magnitude lower than during the summer (e.g. [START_REF] Derwent | Seasonal cycles in shortlived hydrocarbons in baseline air masses arriving at Mace Head, Ireland[END_REF], the dominant NO X sink is conversion of dinitrogen pentoxide (N 2 O 5 ) to HNO 3 on aerosol. This occurs when photolysis is low, allowing NO 3 (formed from the reaction of NO 2 with O 3 ) to build up. NO 3 reacts with NO 2 to form N 2 O 5 . The reaction of N 2 O 5 with H 2 O is slow in the gas-phase [START_REF] Tuazon | The reaction of gasphase N 2 O 5 with water vapour[END_REF] but occurs rapidly in aerosol.

A modelling study by [START_REF] Dentener | Reaction of N 2 O 5 on Tropospheric Aerosols: Impact on the Global Distributions of NO x , O 3 , and OH[END_REF] predicted that changes to the loss of NO x via sulfate aerosol could have a significant effect on northern hemisphere NO x concentrations and that these changes would also affect O 3 and OH concentrations. Subsequent modelling studies have broadly agreed with the magnitude of the changes suggested by Dentener and

Crutzen [START_REF] Brown | Nighttime radical observations and chemistry[END_REF]. Based on these studies, the large decrease in sulfate aerosol observed in Europe and the United States since 1980 [START_REF] Berglen | Sulphate trends in Europe: Are we able to model the recent observed decrease?[END_REF][START_REF] Turnock | Modelled and observed changes in aerosols and surface solar radiation over Europe between 1960 and 2009[END_REF] (driven by a 60 -70 % decrease in SO 2 emissions [START_REF] Smith | Anthropogenic sulfur dioxide emissions: 1850 -2005[END_REF]) would be expected to have led to a decrease in NO X removal by N 2 O 5 hydrolysis, and hence to an increase in the NO X lifetime and atmospheric [NO X ]. This time period is broadly in line with the derived steep increase in the [NO]/[HO 2 ] ratio.

The Hydroperoxyl Radical, HO 2

The complex nature of the interconnectedness of the chemistry make it difficult to predict how HO 2 concentrations will have changed in response to changes in primary emissions.

Changes to Photochemical Oxidation

An alternative explanation for the observed alkyl nitrate trends is that the amount of photochemical processing of the air mass changed during the time period. An increase in processing could be caused by a change in either the concentration of the OH radical (assuming photolysis to have remained constant), or by an increase in the transport time of the air mass from the source region to the Arctic. (E4)

Figure 6 The trend in [OH]t calculated using Equation E4 for each of three alkane-alkyl nitrate pairs. This assumes that the ratio [NO]/[HO 2 ] remained constant at 0.5.

Air mass transport time to the Arctic

The transport time, t, of pollutants to the Arctic from source is dependent on (i) the atmospheric transport patterns, and (ii) the source regions of the pollutants. Concerning (i), [START_REF] Kahl | Decadal-scale variability in long range atmospheric transport to the Summit of the Greenland Ice Sheet[END_REF] have suggested that there is a decadal scale (4 -14 years) variability in transport patterns of pollutants from the NH to the Arctic but note no long term trend. [START_REF] Hirdman | Long-term trends of black carbon and sulphate aerosol in the Arctic: changes in atmospheric transport and source region emissions[END_REF] note that while changes to transport patterns can explain much of the inter-annual variability of Arctic concentrations of black carbon and sulfate aerosol (pollutants with similar source regions to the alkanes), they played only a minor role in long term changes. [START_REF] Eckhardt | The North Atlantic Oscillation controls air pollution transport to the Arctic[END_REF] have

shown that transport of pollutants to the Arctic from European and US source regions is more rapid during positive phases of the weather pattern, the North Atlantic Oscillation (NAO).

During the period 1960 -1980 the NAO was predominantly in a negative phase in winter, between 1980 and 2000 it was predominantly in a positive phase, and since 2000 neither phase has been dominant [START_REF] Hurrell | North Atlantic climate variability: The role of the North Atlantic Oscillation[END_REF]. This suggests that there was more rapid transport of pollutants to the Arctic during the period 1980 -2000 compared to the preceding and succeeding periods. This would mean a shorter processing time for the air masses and hence less alkyl nitrate production and lower alkyl nitrate to alkane ratios. This is the opposite to what we observe in the firn records, suggesting that changes to the NAO are unlikely to be responsible for the observed alkyl nitrate trends. Concerning (ii), the main sources of anthropogenic emissions to the Arctic of gases with lifetimes on the order of a few weeks, particularly during the winter, have been identified as being northern Eurasia (e.g. [START_REF] Shindell | A multi-model assessment of pollution transport to the Arctic[END_REF][START_REF] Stohl | Arctic smoke -Record high air pollution levels in the European Arctic due to agricultural fires in Eastern Europe in spring 2006[END_REF][START_REF] Klonecki | Seasonal changes in the transport of pollutants into the Arctic troposphere -model study[END_REF]. Emissions from Europe and North America have likely followed a similar declining trend in recent years (e.g. [START_REF] Warneke | Multiyear trends in volatile organic compounds in Los Angeles[END_REF], thus the relative contribution from each region will not have changed dramatically.

The Hydroxyl Radical, OH

The alternative explanation for an increase in photochemical processing is an increase in [OH].

This would represent a regional trend in [OH] representative of only regions from and through which air masses are transported to the Arctic, and would relate primarily to the winter (since the alkyl nitrate and alkane signals in the firn are dominated by winter time concentrations). It is noted that the increased chemical processing observed could also result from an increase in an oxidant other than OH, e.g. atomic chlorine, as suggested in Helmig et al. (2014b). Studies using changes to atmospheric mixing ratios of methyl chloroform (CH 3 CCl 3 ) have suggested that global mean OH concentrations are 'well buffered' (e.g. [START_REF] Montzka | Small Interannual Variability of Global Atmospheric Hydroxyl[END_REF].

Since the main sink of CH 3 CCl 3 in the atmosphere is reaction with OH, and the emission sources and other sinks are thought to be well constrained, the variation in its observed mixing ratios at a number of remote sites can be used to infer variations in global [OH]. Global mean

[OH] has been inferred in this way in a number of studies [START_REF] Prinn | Atmospheric trends and lifetime of CH 3 CCl 3 and global OH concentrations[END_REF][START_REF] Prinn | Evidence for substantial variations of atmospheric hydroxyl radicals in the past two decades[END_REF][START_REF] Prinn | Evidence for variability of atmospheric hydroxyl radicals over the past quarter century[END_REF][START_REF] Rigby | Renewed growth of atmospheric methane[END_REF][START_REF] Montzka | Small Interannual Variability of Global Atmospheric Hydroxyl[END_REF]. The most recent of these [START_REF] Montzka | Small Interannual Variability of Global Atmospheric Hydroxyl[END_REF] and 1991 but this increase has been questioned in more recent work (e.g. [START_REF] Montzka | Small Interannual Variability of Global Atmospheric Hydroxyl[END_REF].

However, there are a growing number of observational data sets of trace gases in the NH which show trends since 1980 that could be explained, at least in part, by changes to the concentration of the OH sink. E.g. decreasing Arctic alkane mixing ratios (Helmig et al., 2014b;Aydin et al., 2011); decreasing Arctic CO mixing ratios [START_REF] Petrenko | A 60 yr record of atmospheric carbon monoxide reconstructed from Greenland firn air[END_REF]; increasing d 13 C of methane [START_REF] Monteil | Interpreting methane variations in the past two decades using measurements of CH 4 mixing ratio and isotopic composition[END_REF][START_REF] Sapart | Can the carbon isotopic composition of methane be reconstructed from multisite firn air measurements?[END_REF]; decreasing dC 16 O of Arctic CO [START_REF] Wang | The isotopic record of Northern Hemisphere atmospheric carbon monoxide since 1950: Implications for the CO budget[END_REF].

A recent multi-model inter-comparison exercise of seventeen global chemical transport models, showed agreement for a small increasing trend in global mean [OH] of 3.5 (± 2.2) % between 1980 and 2000 and a slightly larger [OH] increase in the northern hemisphere of 4.6 (± 1.9) % [START_REF] Naik | Preindustrial to present-day changes in tropospheric hydroxyl radical and methane lifetime from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP)[END_REF]. Dalsøren et al. (2015) determined an increase in global mean [OH] of about 10 % between 1970 and 2006 from modelled increases of the methane lifetime.

Discussion

The alkyl nitrate trends presented herein suggest a profound change to the chemical state of the northern hemisphere mid-high latitude atmosphere in winter between the 1970s and the late 1990s and then again between the late 1990s and the mid-2000s.

A key species of the tropospheric chemistry cycle, tightly linked to NO X and HO X , is ozone (Figure 1). Ozone mixing ratios increased at sites across the NH during the second half of the twentieth century, roughly doubling since 1960 [START_REF] Parrish | Long-term changes in lower tropospheric baseline ozone concentrations at northern mid-latitudes[END_REF]. Ozone production is This work also implies that there may have been a change in [OH]. Indeed due to the connected nature of the chemistry of NO X , ozone and OH it seems unlikely, given the implied increases in NO X suggested here, and the recorded increases in ozone, that there was not a commensurate increase in OH during this period. Furthermore, the primary production of OH, denoted by P(OH), is dependent on ozone mixing ratios, either via photolysis of ozone (Equation 5- [START_REF] Smith | Concentrations of OH and HO 2 radicals during NAMBLEX: Measurements and steady state analysis[END_REF] or via ozonolysis of alkenes [START_REF] Johnson | The gas-phase ozonolysis of unsaturated volatile organic compounds in the troposphere[END_REF].

P(OH) = 2f[O 3 ] x j(O 1 D) (E5)
Where f is the fraction of O( 1 D) that reacts with water vapour. Ozone has increased at background sites between 1960 and 2000 [START_REF] Parrish | Long-term changes in lower tropospheric baseline ozone concentrations at northern mid-latitudes[END_REF] and measured water vapour has also increased slightly [START_REF] Hartmann | Observations: Atmosphere and Surface[END_REF]. It therefore seems that the primary production of OH in the background atmosphere must have increased through the final decades of the past century. On the other hand, the primary sink of OH in the background atmosphere, CO, has decreased by about 15% since 1980 [START_REF] Petrenko | A 60 yr record of atmospheric carbon monoxide reconstructed from Greenland firn air[END_REF], with the secondary sink, CH 4 , having increased between 1980 and 2000 by about 15%.

The global growth rate of methane in the atmosphere continually declined throughout the period of the 1970s to 2000, culminating in the 'methane pause ' between 1999' between and 2006' between (Dlugockencky et al. 2009)). A possible cause of this change in growth rate in methane could be an increase in OH concentration (e.g. [START_REF] Dalsøren | Atmospheric methane evolution the last 40 years[END_REF]. This is consistent with an increase in [OH] also being the cause of the trend in the ratio of alkyl nitrate to parent alkane seen in this work. It should be noted, though, that the majority of the OH oxidation of methane occurs in the tropics (e.g. [START_REF] Bloss | The oxidative capacity of the troposphere: Coupling of field measurements of OH and a global chemistry transport model[END_REF], while any increase in OH suggested by the work herein must be viewed as representative only of the mid-high latitude northern hemisphere and the winter time. However, the processes suggested herein, such as changes to the N 2 O 5 sink, have been shown to be effective at a hemispheric scale. In a recent inter-model comparison project, ACCMIP, it was shown that models failed to capture the measured magnitude of the increase in ozone over recent decades [START_REF] Parrish | Long-term changes in lower tropospheric baseline ozone concentrations: Comparing chemistry-climate models and observations at northern mid-latitudes[END_REF], in particular the steep increase seen between 1980 and 2000. This failure to capture measured changes to ozone may suggest that models are likely to under-estimate changes to OH production, from ozone photolysis or reactions of ozone with alkenes, over the same period (i.e. they may be larger than the 4.6 (± 1.9) % reported in [START_REF] Naik | Preindustrial to present-day changes in tropospheric hydroxyl radical and methane lifetime from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP)[END_REF] for the NH).

Including alkyl nitrate chemistry and using the alkyl nitrate measurements presented herein could provide a valuable constraint for modelling changes to NOx and HOx back to the middle of the twentieth century in global chemical transport models.

Figure 1

 1 Figure 1 Schematic of the NO x -O 3 -OH relationship in the background troposphere.

  Atmos.Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.

  model input. It is noted that the latter part (post-1995) of the model derived scenarios for 2+3pentyl nitrate is rather sensitive to the inclusion or exclusion of the measurement at 34.72 m (the most shallow measurement used). The scenarios presented in this work are based on including this measurement. Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.

Figure 2

 2 Figure 2 also shows the atmospheric histories of the same three alkyl nitrates derived from firn air fromNorth GRIP, Greenland, up to 2001, presented in Worton et al. (2010) (pink shaded area). There is very good agreement between the derived trends at the two sites. Differences can be attributed to the limited number of measurements at both sites, possible drift in the calibration standard used and uncertainties in the firn modelling. Both sites show the same large increase in mixing ratios from the 1970s to the 1990s. Importantly, the NEEM records show that the turnover and subsequent decline in mixing ratios, the beginnings of which were evident in the North GRIP records, appears to continue through the 2000s. However, as noted in Section 2.3, the derived atmospheric history of 2+3-pentyl nitrate is sensitive to the inclusion of the measurement at 34.72 m. The scenario that did not include this measurement was almost flat from 1995 to 2008 rather than declining as in Figure2.There are very limited in-situ measurements of alkyl nitrates in the Arctic and even fewer that cover a whole seasonal cycle.[START_REF] Swanson | Seasonal variation of C2-C4 nonmethane hydrocarbons and C1-C4 alkyl nitrates at the Summit research station in Greenland[END_REF] report the seasonal cycle of 2-butyl nitrate at the Summit station,Greenland (72.34 N, 38.29 W, 3250 m a.s.l), from June 1997 to June 1998 based on samples collected roughly every two days.[START_REF] Dibb | An overview of air-snow exchange at Summit, Greenland: Recent experiments and findings[END_REF] report monthly mean measurements of 2-butyl nitrate for the period June 2000 through to August 2002 based on samples taken roughly weekly also from Summit. In order to compare these in-situ

Figure 2

 2 Figure 2 Concentration-depth profiles in the firn and the model derived atmospheric histories. Right panel: The concentration-depth profiles measured in the firn (ppt): red open circles: measured mixing ratios (ppt) with error bars indicating the 2-σ uncertainty; solid blue line: best fit of the firn model, dashed blue lines indicate the 2-σ combined analytical and model uncertainties. Left panel: Atmospheric histories of the alkyl nitrates derived from the firn air measurements using the inverse modelling technique described within (solid blue lines). Dashed lines represent the 2-σ confidence margins of the model calculations, combining the analytical and model uncertainties. Pink shaded area: atmospheric histories presented in Worton et al. (2012) derived from firn air measurements at North GRIP, Greenland. Red open circles: Annual average of in-situ measurements at Summit, Greenland (see text for details) with 1-σ uncertainty.

Figure 3

 3 Figure 3 Atmospheric histories of 2-butyl nitrate (red) and its parent alkane, n-butane (blue), derived from firn measurements at NEEM, Greenland.

4Figure 4

 4 Figure4The trend in the [NO]/[HO 2 ] ratio of the air masses in which the alkyl nitrates were formed, calculated using Equation E3 for each of three alkyl nitrate/alkane pairs. This assumes that the amount of photochemical processing, [OH]t, remained constant at 5 × 10 11 molecules cm -3 s.

Figure 4

 4 Figure 4 shows the trend in the [NO]/[HO 2 ] ratio derived using Equation E3. A value of 5 × 10 11 molecules cm -3 s is used for the constant [OH].t. This is based on a mean transport time of air masses from Europe (from where the majority of winter-time pollutants are transported to the Arctic -see Section 6.1) to the Arctic in the winter of ten days (Stohl, 2006), and a mean winter-time [OH] of ~ 6 × 10 5 cm -3 (in reasonable agreement with that derived by Derwent et al. (2012) for the North Atlantic in winter-time). Halving the assumed [OH].t to 2.5 × 10 11 molecules cm -3 s increases the magnitude of the relative changes by a factor of ~ 2. Whereas doubling [OH].t, only decreases the derived changes by ~ 15 %.

Figure 5

 5 Figure5shows how NO X emissions from North America, Europe and Russia have varied between 1970 and 2008, taken from the bottom-up estimates of the EDGAR database (EDGAR v4.2, http://edgar.jrc.ec.europa.eu). Emissions were fairly constant between 1970 and 1990 and then fell by about 25 % from 1990 to 2008. This is in good agreement with the NO X emission trends for OECD Europe presented in[START_REF] Vestreng | Evolution of NO x emissions in Europe with focus on road transport control measures[END_REF].

Figure 5

 5 Figure 5 Blue solid line (left axis): The trend in NO X emissions (Tg yr -1 ) from North America, Europe and Russia for the period 1970 to 2008 (EDGAR v4.2, http://edgar.jrc.ec.europa.eu). Green points (and dashed line) (right axis): OECD Europe NO X emissions (Tg yr -1 ) from Vestreng et al. (2009).

  Atmos.Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.

  photochemical processing, [OH]t, can be calculated using the changes to the [RONO 2 ]/[RH] ratio (assuming that the [NO]/[HO 2 ] ratio has remained constant).

Figure 6

 6 Figure 6 shows the trends in [OH].t derived from the alkyl nitrate-alkane pairs if a constant value for [NO]/[HO 2 ] is assumed. The value used for the constant [NO]/[HO 2 ] was 0.5. This is in line with the [NO]/[HO 2 ] derived using Equation E3 with a fixed [OH]t value of 5 × 10 11 molecules cm -3 s.

Figure 6

 6 Figure 6 shows that the observed [RONO 2 ]/[RH] ratios between around 1970 and the late 1990s could be explained by a relative change in [OH]t of a factor of between 2.4 (2-butyl nitrate) and 3.3 (3-methyl-2-butyl nitrate). Increasing the [NO]/[HO 2 ] value by an order of magnitude, to 5, decreases the derived [OH]t changes by ~ 25 %. Decreasing the [NO]/[HO 2 ] value by an order of magnitude, to 0.05, increases the derived [OH]t changes by ~ 5 %.

  Atmos.Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.

  Atmos.Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License. positively linked to [NO X ] (at low NO X concentrations such as the background atmosphere).Hence, an increase in the [NO]/[HO 2 ] ratio from around 1970 to the mid-1990s is consistent with long-term trends in ozone in the background atmosphere. Furthermore, while the alkyl nitrate measurements represent changes to the winter-time atmosphere, the ozone trends are seen in both summer and winter. If these are being driven by increases to [NO X ] in the background atmosphere, then this suggests that the chemical changes to the atmosphere may have been present throughout the year and are not just a winter time phenomenon.

  Atmos.Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.

  as those presented here are fundamental to improving our understanding of trends in atmospheric composition during the twentieth century. The long-term trends of alkyl nitrates presented herein suggest major changes to the chemical state of the atmosphere during the past five decades. The trends were likely driven by changes to the [NO]/[HO 2 ] ratio during this period. A 4 -5 fold increase in the [NO]/[HO 2 ] ratio is suggested between the 1970s and late 1990s. This is not consistent with reported changes to NO X emissions, but may have been driven by a reduction in the NO X sink. The recent decreases (since the late 1990s) in alkyl nitrate concentrations are in qualitative agreement with recent decreases in NO X emissions and in measured NO X concentrations. Changes to concentrations of the hydroxyl radical, or to the transport time of air masses from source regions to the Arctic may also have contributed to the observed trends. of Polar Programs) and the EU Seventh Framework programs. We are indebted to Jakob Schwander of the Physics Institute at the University of Bern, Switzerland for collecting the firn air samples at NEEM, and Thomas Blunier of the Centre for Ice and Climate at University of Copenhagen, Denmark for leading the NEEM gas consortium. We thank Chelsea Thompson for useful discussions.

  This is a minor channel of the RO 2 + NO reaction (Reaction R14a) which generally leads to ozone production via recycling of NO to NO 2 and the subsequent photolysis of NO 2 . RO 2 can also react with HO 2 (the hydroperoxyl radical) (Reaction R11) to

	deposition.
	RH + OH à RO 2	(R13)
	RO 2 + NO à RCHO + NO 2	(R14a)
	à RONO 2	(R14b)
	RO 2 + HO 2 à ROOH	(R11)
	RONO 2 + OH à products	(R15)
	RONO 2 à products	(R16)
	2 Methodologies
	2.1 Firn Sampling
	Firn air samples were collected at the NEEM site, Greenland (77.45°N, 51.07°W, 2484m a.s.l)
	from two boreholes between 14 th and 30

form a peroxide (ROOH). The probability of RO 2 reacting with NO (leading to alkyl nitrate production) is thus governed by the ratio [NO]/[HO 2 ]. Alkyl nitrates are lost from the atmosphere by reaction with OH (Reaction R15), photolysis (Reaction R16) and wet/dry th July 2008 ("2008 EU hole" and "2008 US hole"). Further samples were collected from another NEEM borehole during July 2009 ("2009 hole").

  reported little inter-annual variability in mean global atmospheric [OH] estimating roughly 5% variation from the mean value during the period 1997 -2008, but this does not cover the period of interest here (1970 to later 1990s). Similarly, Kai et al. (2011) inferred a low variability in global [OH]between 1998 and 2005 based on a relatively constant δ-D-CH 4 . Earlier work using methyl chloroform[START_REF] Prinn | Evidence for substantial variations of atmospheric hydroxyl radicals in the past two decades[END_REF] reported an increase in NH[OH] of roughly 40% between 1979

Acknowledgements

This work was supported by funding from the UK Natural Environment Research Council (NE/F021194/1 & NE/M003248/1). NEEM is directed and organized by the Centre of Ice and Climate at the Niels Bohr Institute and US NSF, Office of Polar Programs. It is supported by funding agencies and institutions in Belgium (FNRS-CFB and FWO), Canada (NRCan/GSC), China (CAS), Denmark (FIST), France (IPEV, CNRS/INSU, CEA and ANR), Germany (AWI), Iceland (RannIs), Japan (NIPR), South Korea (KOPRI), The Netherlands (NWO/ALW), Sweden (VR), Switzerland (SNF), the United Kingdom (NERC) and the USA (USNSF, Office Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-927, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 30 November 2016 c Author(s) 2016. CC-BY 3.0 License.