Secular change in atmospheric Ar/N 2 and its implications for ocean heat uptake and Brewer-Dobson circulation

. Systematic measurements of the atmospheric Ar/N 2 ratio have been made at ground-based stations in Japan and Antarctica since 2012. Clear seasonal cycles of the Ar/N 2 ratio with summertime maxima were found at middle to high latitude stations, with seasonal amplitudes increasing with increasing latitude. Eight years of the observed Ar/N 2 ratio at Tsukuba (TKB) and Hateruma (HAT), Japan showed interannual variations in phase with the observed variations in the global ocean heat content (OHC). We calculated secularly increasing trends of 0.75±0.30 and 0.89±0.60 per meg yr -1 from the Ar/N 2 ratio 15 observed at TKB and HAT, respectively, although these trend values are influenced by large interannual variations. Sensitivity test by using a 2-dimensional model with the Brewer-Dobson circulation (BDC) scenarios indicated the possibility of the secular trend in the surface Ar/N 2 ratio being modified significantly by the gravitational separation in the stratosphere. The secular trend of the Ar/N 2 ratio at TKB, corrected for gravitational separation under the assumption of weakening of BDC simulated by the 2D model, was 0.60±0.30 per meg yr -1 . By using a conversion factor of 3.5x10 -23 per meg J -1 by assuming a 20 1-box ocean with a temperature of 3.5 °C, an average OHC stations examined the relationship of Ar and N 2 fluxes with air-sea heat flux. prepared gravimetric standard mixtures of Ar, O 2 , N 2 and CO 2 . SS, YT, DG, KI, YN, NA and SM examined the results and provided feedback on the manuscript. All the authors 440 approved the final manuscript.


Responses to Referee 2:
Thank you very much for your significant and useful comments on the paper "Secular change in atmospheric Ar/N2 and its implications for ocean heat uptake and Brewer-Dobson circulation" by Ishidoya et al. We have revised the manuscript, considering your comments and suggestions. Details of our revision are as follows; The manuscript is challenging to review as it touches on a wide range of topics but information provided is often (too) sparse. As far as I can see, the main points of the manuscript are: (i) Ar:N2 measurements are accurate enough to detect trends, which was not possible previously (Line 37).
(ii) Variations in the ratio over a decade are well correlated with OHC variations as reported by NOAA/NCEI.

(iii) There is a quantitative mis-match if only OHC is considered.
(iv) A model simulation with imposed change in the stratospheric BD circulation suggests that the near surface Ar:N2 ratio is as sensitivite to the BD circulation as to OHC variations (for the observed magnitude of changes in either of these).
(v) If one assumes that the OHC measurements and their impact on Ar:N2 are correct, one infers that over the period of measurements (2012-2019) the BD circulation slowed down.
The authors discuss the last point in some detail, and seem to be concerned that the sign of the change in the BD circulation is the opposite of a supposedly long term strengthening of the BD circulation. There is no reason, however, why a short period such as the one studied here should show the long term trend given the large interan nual variability of the stratospheric circulation. I recommend to shorten this somewhat misleading discussion (the authors are most likely not observing "trends"), and instead strengthen the points (i)-(iv).
Thank you for your comments and suggestion. As you pointed out, the observational period of 8-years in the present study is not enough to discuss the long-term trend of the BDC. In addition, the troposphere, and not only the ocean, does not mix perfectly on a timescale of a year, and that the surface-ocean temperature anomalies would be a large source of interannual variation on a yearly timescale for the observed d(Ar/N2) in the near-surface air. Therefore, we recognize the possibility that the secular trends observed at the surface and simulated by the 2-D model do not represent a long-term trend but is a part of the large interannual variation. Nevertheless, it would of interest to see if it is possible to obtain a scientifically "meaningful" OHC change based on the firstly-reported secular δ(Ar/N2) trend at the surface and the new concept of δ(Ar/N2)Ω. Taking these into consideration, we have revised the related sentences throughout the manuscript, and stated clearly in the revised sentences that the 2-D model simulations are not enough to represent the mechanism to drive AoA change in the real atmosphere but it does serve as a kind of sensitivity test (the revised sentences have been highlighted by blue color : for example, lines [216][217][218][219][220][221][222][223][224][225][226][227][228][229][230][231][232][293][294][295][296][297][298][299][300][301][302][303][304]and Appendix A and B). Thank you for your comments. As you expected, another reviewer commented on this matter and we have added the sentences not only to describe the uncertainties for seasonal cycles (lines 154-161) but also to show the detail method to extract the secular trend of δ(Ar/N2) observed at the surface stations (lines 216-232).

(b) The time filtering and frequency separation method needs to be better explained -based on the information provided in this manuscript, it is not possible to
reconstruct what precisely the authors have done. This is an extremely important part of the argument as the claim is that seasonal variations are due to local SST variations, whereas the longer-term changes are due to global OHC changes (i.e. the seasonal cycle at TKB and HAT differ (Fig 3), but their longerterm variations (Fig 4) are highly correlated). Also, the analysis then uses largely only TKB and HAT -but why not also show COI in Figure 4? Please use 2 different colors for AHT and TKB in Figure 4a -the separation with "+" and "o" does not allow to see the differences (Looking at Figure 2; it is surprising that these two timeseries align as well as shown in Figure 4).
Lines 140-147 and Fig. 4: The sentences have been added to explain the time filtering and frequency separation method more in detail. Also, we have added the data at COI in A1 and A2). We have also added the sentences and figure to discuss the annual mean meridional distribution of the AoA trend (yrs yr -1 ) calculated using the SOCRATES model for weakened BDC simulation (lines 421-430 and Fig. A3). We have recognized the simulations using the 2-D model in the present study, made by arbitrarily changing the MMC only with fixed horizontal mixing, are not enough to represent the mechanism to drive AoA change in the real atmosphere but it does serve as a kind of sensitivity test. Therefore, we have clarified the limitation of the 2-D model (lines 293-304). We consider our simulation is worthwhile as the sensitivity test since it constitutes the first step to investigate the effect of gravitational separation of the whole atmosphere on the surface Additional points 1) Figure 5: My understanding is that this is simply an idealized calculation for 10 years; if so, please change x-label (labels 2000 -2010 suggest that this is specifically for this period; which is confusing also because OHC changes are taken from this period, see next item).

4) Caption for Fig. 2
We have removed the incorrect sentence "All data are corrected for the scale drift of the primary standard air shown in Fig. 1 (b)." in the caption for Fig. 2 in ACPD paper since we have not applied such correction to the data.

Responses to Referee 3:
Thank you very much for your significant and useful comments on the paper "Secular change in atmospheric Ar/N2 and its implications for ocean heat uptake and Brewer-Dobson circulation" by Ishidoya et al. We have revised the manuscript, considering your comments and suggestions. Details of our revision are as follows; The one minor comment I have is that the one-box ocean model is somewhat overemphasized, and perhaps the description could be simplified and shortened, because we know very well that the troposphere does not perfectly mix on a timescale of a year, nor does the ocean mix perfectly on a timescale of one year, so it is not really surprising that a one-box ocean model fails to match the observations of Ar/N2 at surface stations around Japan. I understand that the authors constructed the onebox model as a "straw man", to be shot down, but they could greatly simplify and shorten the discussion, while making it clear that they do not expect to be able to match their observations with a one-box model of the ocean. In fact, some readers may be confused, the way the authors have written about this one-box model (they seem to imply that they expected it to be able to match their observations). But of course the surface-ocean temperature anomalies are by far the largest source of noise on yearly timescales, for tracers like Ar/N2, measured in near-surface air.
They should clarify that a true tropospheric-average Ar/N2, measured from aircraft (which is of course too expensive and so is prohibitive), would be needed to actually compare Ar/N2 observations with modeled Ar/N2. (So in some sense they are raising a false expectation with their one-box model exercise.) On the whole, this is a groundbreaking and important paper, and will make a fine contribution to ACP. This is clearly excellent, highest-quality work! Lines 25-28, 196 and 205-215: Thank you very much for giving high evaluation to the present study. As you pointed out, the one-box ocean model is not enough to evaluate responses of the atmospheric δ(Ar/N2) to changes in the air-sea heat flux in detail. Therefore, we have stated the limitation of the one-box ocean model clearer. We have also added a reference showing the renewal time of permanent pycnocline water in the North Pacific, which would make the readers imagine that the secular trend of the δ(Ar/N2) for 8-years in the present study mainly reflects the OHC change except deep ocean. We recognize the OHC changes estimated in the present study, based on the observed 8-years δ(Ar/N2) trend combined with the 2-D atmospheric model and the onebox ocean model, are insufficient to suggest some revisions of the OHC from ocean temperature measurements. Nevertheless, it would of interest to see if it is possible to obtain a scientifically "meaningful" OHC change based on the firstly-reported secular δ(Ar/N2) trend and the new concept of δ(Ar/N2)Ω.
Line 39: The words "Argo float" have been changed to "Argo floats", as suggested.
(2) line 47 -". . .since it has been believed that the gravitational separation near the surface is too small to be detected..." This isn't really accurate, in the sense that there is no need for the gravitational settling process IN THE troposphere to be significant, in order for the stratosphere to affect the troposphere. Perhaps instead you could say something like, ". . .since it has been believed that the gravitational settling in the stratosphere is fairly small and constant in time, along with the fact that the troposphere has 10x more molecules than the stratosphere." Or perhaps you meant to say, ". . .the gravitational separation signal from the stratosphere is too small to be detected at the surface."? This is accurate.
Line 50-51: The words since it has been believed that the gravitational separation near the surface is too small to be detected" have been changed to "since it has been believed that the gravitational separation signal from the stratosphere is too small to be detected at the surface", as suggested.
(3) line 51 ". . .long-term changes in the Ar/N2 ratio near the surface are expected to be extremely small.." Line 54: The word "near the surface is" have been changed to "near the surface are". Line 62: The words "trend of Ar/N2 ratio" have been changed to "trend of the Ar/N2 ratio" (6) line 60 "Atmospheric Ar/N2 has been observed. . ..." (there is no need to include "ratio" here) Line 65: We removed the word "ratio", as suggested. Lines 77-79: The sentence has been modified, considering your suggestion. Line 105: The words "glass flask" have been changed to "glass flasks".
(9) Line 102 "per meg units as follows." Line 106: The words "per meg unit" have been changed to "per meg units". Lines 115-116: The words "the uncertainty of" have been changed to "an uncertainty of".
(12) Line 122 "..but they did not correlate.." Line 125: The words "but did not" have been changed to "but they did not".
Line 124 ". . . must have been superimposed. . ." Line 127: The words "must have superimposed on" have been changed to "must have been superimposed on".
(14) Line 126 ". . . due to a temperature. . ." Line 129: The words "due to temperature" have been changed to "due to a temperature".
(15) Line 138 ". . . by a fundamental sine-cosine. . ." Line 141: The word "fundamental" have been changed to "a fundamental sine-cosine". Line 151: The word "enhancing" have been changed to "due to". Line 162: The word "increase" has been changed to "increases". Line 166: The words "found seasonal" have been changed to "found a seasonal". Line 170: The words "than 14±6" have been changed to "than the 14±6". The OHC values are shown as anomalies from the baseline value observed in mid-1980s.
In the figure we also plot a interannual variation of the OHC values obtained by using the same digital filtering technique used in Fig. 2, and globally averaged surface temperature anomalies (Japan Meteorological Agency, http://www.data.jma.go.jp/cpdinfo/temp/nov_wld.html).".
Perhaps say "As a first approximation we modeled . . ." Line 197: The words "We boldly modeled" have been changed to "As a first approximation we modeled", as suggested.

(23) line 196 "As mention in the Introduction, "
Line 235: The words "in Introduction" have been changed to "in the Introduction". Lines 239-241: The sentence has been modified considering your suggestion as "Therefore, we need to explore the possibility of tropospheric d(Ar/N2) variations caused by changes in the stratospheric gravitational separation that influence whole troposphere.".
Line 306: The word "seesaw" has been changed to "inverse".
(26) line 265 "inputted" is an awkward word. Perhaps instead use "heat is added to a . . ." Line 310: The words "is inputted" have been changed to "is added to". Line 313: The words "are non-negligible" have been changed to "is a non-negligible trend".
(30) Caption for Fig. 2 We have removed the incorrect sentence "All data are corrected for the scale drift of the primary standard air shown in Fig. 1 (b)." in the caption for Fig. 2 in ACPD paper since we have not applied such correction to the data. Abstract. Systematic measurements of the atmospheric Ar/N2 ratio have been made at ground-based stations in Japan and Antarctica since 2012. Clear seasonal cycles of the Ar/N2 ratio with summertime maxima were found at middle to high latitude stations, with seasonal amplitudes increasing with increasing latitude. Eight years of the observed Ar/N2 ratio at Tsukuba (TKB) and Hateruma (HAT), Japan showed interannual variations in phase with the observed variations in the global ocean heat content (OHC). We calculated secularly increasing trends of 0.75±0.30 and 0.89±0.60 per meg yr -1 from the Ar/N2 ratio 15 observed at TKB and HAT, respectively, although these trend values are influenced by large interannual variations. Sensitivity test by using a 2-dimensional model with the Brewer-Dobson circulation (BDC) scenarios indicated the possibility of the secular trend in the surface Ar/N2 ratio being modified significantly by the gravitational separation in the stratosphere. The secular trend of the Ar/N2 ratio at TKB, corrected for gravitational separation under the assumption of weakening of BDC simulated by the 2D model, was 0.60±0.30 per meg yr -1 . By using a conversion factor of 3.5x10 -23 per meg J -1 by assuming a 20 1-box ocean with a temperature of 3.5 °C, an average OHC increase rate of 17.1±8.6 ZJ yr -1 for the period 2012 -2019 was estimated from the corrected secular trend of the Ar/N2 ratio. This value is consistent with 12.2±1.2 ZJ yr -1 reported by ocean temperature measurements. On the other hand, the OHC increase rate of 25.1±8.6 ZJ yr -1 for the corresponding period, estimated from the secular trend of the Ar/N2 ratio corrected under the assumption of enhanced BDC, is significantly larger than the OHC value from the ocean temperature measurements. Although the effect of the actual atmospheric circulation on 25 the Ar/N2 ratio is still unclear and longer-term observations are needed to reduce uncertainty of the secular trend of the surface Ar/N2 ratio, the analytical results obtained in the present study imply that the surface Ar/N2 ratio is an important tracer for detecting spatiotemporally-integrated changes in OHC and BDC.

Introduction 30
The Ar/N2 ratio of air is a unique tracer for detecting changes in the spatiotemporally-integrated air-sea heat flux or ocean heat content (OHC). This is because variations in the Ar/N2 ratio at the Earth's surface are driven by air-sea Ar and N2 fluxes that reflect changes in the solubility of these gases in seawater (e.g. Keeling et al., 2004). The increase in the global OHC is one of the most important parameters for evaluating earth's climate system (e.g. Trenberth and Fasullo, 2010). The relative temperature dependence of the solubility of Ar is larger than that of N2, so that the atmospheric Ar/N2 ratio increases with 35 increasing ocean temperature. Other noble gases behave in a similar way. Bereiter et al. (2018), for example, analyzed Kr/N2 and Xe/N2 ratios in air trapped in ice cores and estimated that the mean global ocean temperature increased by 2.57±0.24 °C over the last glacial transition (20,000 to 10,000 years ago). In terms of ocean heat content, the present-day global OHC increase is evaluated from analysing the ocean temperature measurements using Argo floats (e.g. Cheng et al., 2019). Therefore, precise measurements of the Ar/N2 ratio in air can be used as an independent validation of the OHC estimation from the ocean 40 data. However, long-term changes in the atmospheric Ar/N2 ratio have never been reported so far due to difficulties in detecting a trend with sufficient accuracy, although a few past studies observed its seasonal variations (Keeling et al., 2004;Blaine, 2005;Cassaer et al., 2008;Ishidoya and Murayama, 2014).
We have reported in past studies that the Ar/N2 ratio decreases with increasing altitude in the stratosphere due to gravitational separation of the atmospheric constituents (e.g. Ishidoya et al., 2013;Sugawara et al., 2018). The magnitude of gravitational 45 separation is determined by a balance between mass-independent atmospheric transport, that is, advection and eddy diffusion, and mass-dependent molecular diffusion in the atmosphere. This implies that gravitational separation will be influenced by the atmospheric circulation changes. Therefore, we can use the observed gravitational separation as an indicator of the Brewer-Dobson circulation (BDC) in the stratosphere. There has been no study to evaluate the effect of gravitational separation changes in the stratosphere on the concentrations and isotopic ratios of atmospheric trace gases at the surface, since it has been believed 50 that the gravitational separation signal from the stratosphere is too small to be detected at the surface. Therefore, in our previous studies, we simulated "d", which is an indicator of gravitational separation derived from the Ar/N2 ratio and stable isotopic ratios of N2, O2 and Ar, using atmospheric transport models by assuming the surface d to be zero (Ishidoya et al., 2013;Sugawara et al., 2018;Belikov et al., 2019). However, because long-term changes in the Ar/N2 ratio near the surface are expected to be extremely small, e.g. 0.00026 % corresponding to a heat input of 100 ZJ (1 zetajoule = 10 21 J) into a 10 °C 55 ocean (Keeling et al., 2004), a very small secular change in the stratospheric gravitational separation signal near the surface may modify the long-term change in the surface Ar/N2 ratio. If so, an evaluation of gravitational separation of the whole atmosphere is needed for a precise estimate of the global OHC increase based on long-term changes in the surface atmospheric Ar/N2 ratio.
In this paper, we present results from an analysis of 8-year-long measurements of the atmospheric Ar/N2 ratio at the ground 60 surface stations and propose "dW" as a new indicator of gravitational separation of the whole atmosphere by using a 2-D model. By using the simulated dW, we derived the secular trend of the Ar/N2 ratio corrected for gravitational separation changes associated with the BDC change. Finally, we estimated the global OHC change based on the corrected Ar/N2 ratio.

Experimental procedures
Atmospheric Ar/N2 has been observed at Tsukuba (TKB; 36°N, 140°E), Japan continuously since February 2012 (Ishidoya 65 and Murayama, 2014). Located on the roof of a laboratory building of the National Institute of Advanced Industrial Science and Technology (AIST), sample air is taken from an air intake by using a diaphragm pump with gas velocity higher than 5 ms -1 (4 mm ID and a flow rate of 4 L min -1 ) at the tip of the air intake to prevent thermally-diffusive inlet fractionation (Sturm et al., 2006;Blaine et al., 2006). The sample air is introduced into a 1 L stainless steel buffer tank after water vapour in the air is reduced by using an electric cooling unit at 2°C. The gas is then exhausted from the buffer tank at a flow rate of about 4 L min -70 1 . A small portion of this exhausted gas is introduced into a 1/8-inch OD stainless-steel tube and any remaining water vapour is removed using a cold trap at -90 °C. Finally, the remaining sample air is vented through an outlet path at a rate of about 10 mL min -1 , and only a miniscule amount of it is transferred to the ion source (or waste line) of a mass spectrometer (Thermo Scientific Delta-V) through an insulated thin fused-silica capillary. As for the reference air, it is always supplied from a highpressure cylinder at a flow rate of about 4 mL min -1 , and a miniscule amount of it is introduced into the ion source (or waste 75 line) of the mass spectrometer through another fused-silica capillary. For the continuous measurements, alternate analyses of the sample and reference air are made repeatedly. The measurement time required to obtain one data value is 62 seconds, but we usually use the 550-data averaged value as the reported Ar/N2 ratio obtained from the continuous observations (about 11hour averaged value). We also measure stable isotopic ratios of N2, O2 and Ar, and O2/N2 ratio and CO2 mole fraction simultaneously, and use the 550-data averaged values for the stable isotopic ratios and one data value without averaging for 80 O2/N2 ratio and CO2 mole fraction; these measurements constitute the continuous observations. Details of the continuous measurement system used are given in Ishidoya and Murayama (2014).
We have also collected air samples at a rate of once per month at Hateruma Island (HAT; 24°N, 124°E) and Cape Ochiishi (COI; 43°N, 146°E), Japan since July 2012 and October 2013, respectively, and at Syowa station (SYO; 69°S, 40°E) since January 2016, for the analyses of the atmospheric Ar/N2 ratio. Each air sample is taken from an air intake by using a diaphragm 85 pump at a flow rate of about 5 L min -1 and filled into a 1 L stainless steel flask whose inner walls are silica-coated after removing water vapour using a cold trap (-40 °C at HAT and COI, and -80 °C at SYO). Similar to our previous study (Ishidoya et al., 2016), the 1 L stainless steel flasks are equipped with two metal-seal valves on each side to equalize the inner pressure to the pressure between the two metal-seal valves to prevent a mass-dependent fractionation due to small leak through the valve. During air sampling, the inner pressure of the flask is kept at an absolute pressure of 0.2 MPa using a backpressure 90 valve (Tohjima et al., 2003). The sample air collected in the flasks are then sent to AIST and analyzed by using the same mass spectrometer as described above. The sample air is supplied from the flask at a flow rate of about 4 mL min -1 through a cold trap (about -50 °C), and a miniscule amount of it is introduced into the ion source (or waste line) of the mass spectrometer through a fused-silica capillary. It is noted that we also used 760 mL glass flasks with a Viton O-ring seal valves on each side for collecting air samples at HAT and COI prior to September, 2015 and January, 2019, respectively. However, we found 95 slight seasonally dependent differences in the Ar/N2 and O2/N2 ratios between the analytical results from the 1 L stainless steel flasks and those from the 760 mL glass flasks. It is interesting to note that the O2/N2 ratio from the 1 L stainless steel flasks agree well with the O2/N2 ratio reported by the National Institute for Environmental Studies (NIES) (e.g. Tohjima et al., 2019), considering the difference in the span sensitivity of the O2/N2 ratio between the AIST and NIES (our unpublished data). Therefore, we have decided to adopt the Ar/N2 data obtained from the 1 L stainless steel flasks and correct the data from the 100 760 mL glass flasks based on the comparison of the Ar/N2 ratios measured from the stainless steel flasks and the glass flasks at HAT for the period October, 2015 -January, 2019. Cause(s) of such an offset between the stainless steel flasks and glass flasks have not been determined yet, but it may be related to the seasonal difference in the ambient temperature during the time the flasks were shipped from the observational sites to our laboratory and its effect on the condition of Viton O-ring seal valves used in the glass flasks. 105 The Ar/N2 ratio is usually reported in per meg units as follows.
Here, the subscripts 'sample' and 'standard' indicate the sample air and the standard gas, respectively. Because Ar constitutes 9,334 μmol mol −1 of air by volume (Aoki et al., 2019), 5 per meg of d(Ar/N2) is equivalent to about 0.05 μmol mol −1 . In this study, d(Ar/N2) of each air sample was determined against our primary standard air (cylinder No. CRC00045) using the mass 110 spectrometer Thermo Scientific Delta-V. Our air standards, which are classified into primary and secondary, are dried ambient air or industrially-purified air-based CO2 standard filled in 48-L high-pressure cylinders. As shown in Fig. 1, variations in the annual average d(Ar/N2) of our 3 secondary standards are within ±0.9 -±2.2 per meg (±1.6 per meg, on average) and nearly stable for 8 years with respect to the primary standard. Therefore, we allowed an uncertainty of ±1.6 per meg for the annual average d(Ar/N2) observed in the present study associated with the stability of the standard air, which corresponds to an 115 uncertainty of ±0.28 per meg yr -1 for the 8-year-long secular trend. We have also prepared high-precision gravimetric standard mixtures of Ar, O2, N2 and CO2, with standard uncertainties for the Ar and O2 molar fractions of 0.6 to 0.8 μmol mol −1 (Aoki et al., 2019). From the measurements of the gravimetrically-prepared standard mixtures by the mass spectrometer, it was confirmed that the span sensitivity of d(Ar/N2) obtained from the mass spectrometer agreed to within 0.2 % of the gravimetric values of the standard mixtures, in the range from -4,500 to +1,800 per meg of d(Ar/N2) against our primary standard air. 120 In the present study, we use 1-week averaged values from the continuous observation at Tsukuba after implementing the following data selection procedure. First, the d(Ar/N2) values with d 15 N higher than 3.0 per meg were excluded from the analyses. As mentioned above, the d(Ar/N2) and d 15 N values are already 11-hour averaged values. We have found these high d 15 N events accompanied by high CO2 events, as well as occasionally with slightly high d(Ar/N2) events, especially in the winter, but they did not correlate with variations in the isotopic ratio of O2 and Ar. Therefore, it appears that some unspecified 125 interferences of mass 29 (possibly 13 C 16 O and/or fragments of hydrocarbons) rather than the molecular-diffusive separation of 15 N 14 N and 14 N 14 N must have been superimposed on the observational results of d 15 N during these events. It is possible that these high CO2 and d 15 N events occur under a stable atmospheric condition in the winter, so that simultaneously observed d(Ar/N2) may also be modified by local effects such as thermally-diffusive separation of Ar and N2 due to a temperature inversion near the surface (Adachi et al., 2006). Therefore, the threshold value of 3.0 per meg was determined to be reasonable, 130 considering that an average d 15 N value of 1.1±1.7 per meg was observed at HAT for the period October 2015 -January 2020, while the observed CO2 mole fractions over the same period were much closer to the values of the background air than those observed at TKB. After the above data selection procedure, 1-week averaged values of d(Ar/N2) were calculated. The measurement uncertainty of the 1-week averaged values of the continuous observation was estimated to be about ±3 per meg as a standard deviation from the best-fit curve shown in Fig. 2 discussed below, while that of the flask air sample measurements 135 was estimated to be about ±7 per meg from repeated analyses of the same air samples.

Atmospheric d(Ar/N2) observed at ground surface stations in Japan and Antarctica
Figure 2 shows atmospheric d(Ar/N2) observed at COI, TKB, HAT and SYO. Best-fit curves to the data and long-term trends obtained using a digital filtering technique (Nakazawa et al., 1997) are also shown. Using this filtering technique, the average 140 seasonal cycle of d(Ar/N2) was modelled by a fundamental sine-cosine and its first harmonic with the respective periods of 12 and 6 months. The residuals obtained by subtracting the average seasonal cycle from the data were interpolated linearly to calculate the daily values of d(Ar/N2), which were then smoothed by a 26th-order Butterworth filter with a cutoff period of 36 months to obtain the interannual variation. The interannual variation was then subtracted from the original data, and the average seasonal cycle was determined again from the residuals. These steps were repeated until we obtained an unchanging interannual 145 variation. In this study, we have defined an average linear increasing/decreasing trend as the "secular trend" (as in Fig. 5 discussed below).
As seen in Fig. 2, we can distinguish clear seasonal d(Ar/N2) cycles at COI, TKB and SYO and some interannual variations at TKB and HAT. Figure 3 shows the detrended values of d(Ar/N2) and average seasonal d(Ar/N2) cycles observed at all 4 sites.
The seasonal maxima were found in the summertime, which is expected since the sea surface temperatures around the 150 observational sites reach seasonal maximum in the summertime, due to the larger relative temperature dependent solubility of Ar compared to that of N2. The peak-to-peak amplitudes of the seasonal d(Ar/N2) cycles were 21±10, 11±4, 5±10 and 32±9 per meg at COI, TKB, HAT and SYO, respectively. The uncertainties for the amplitudes indicate standard deviations of the detrended values from the average seasonal cycle. The uncertainties at COI, HAT and SYO were found to be not only larger than that at TKB, but also ±7 per meg of the uncertainty originated from the repeated analyses of the same flask air sample. 155 This would be due to the fact that the uncertainty of each analysis value of the standard air (about ±5 per meg, black dots in Fig. 1) is superimposed on the uncertainty of each analysis of the flask air sample and continuous measurement. Therefore, a mean squared error expected for the observational data from the flask air sample is about ±9 per meg, which is comparable to the uncertainties in the seasonal amplitudes in Fig. 3. Corrections to the d(Ar/N2) data at HAT and COI prior to September, 2015 and January, 2019 discussed above (corrections based on the comparison of the d(Ar/N2) values from the stainless steel 160 flasks and the glass flasks), could also be contributing to the uncertainties. Similar increases in the seasonal d(Ar/N2) cycle amplitude with increasing latitude were also observed by Keeling et al. (2004) and Cassar et al. (2008). Cassar et al. (2008) also reported on the seasonal d(Ar/N2) cycle at SYO with a peak-to-peak amplitude of 21±8 per meg. The amplitude at SYO found in this study is consistent with that found by Cassar et al. (2008), within the quoted uncertainties. For La Jolla (LJO; 33°N, 117°W), USA located at a similar latitude as TKB, Keeling et al. (2004) and 165 Blaine (2005) found a seasonal d(Ar/N2) cycle with a peak-to-peak amplitude of about 10 per meg. This agrees with the amplitude at TKB observed in this study. Seasonal minima and maxima at SYO and LJO reported by Cassar et al. (2008) and Keeling et al. (2004), respectively, are in general agreement with those observed at SYO and TKB in this study. On the other hand, the seasonal d(Ar/N2) cycle at HAT was not so clear but the average peak-to-peak amplitude may be slightly smaller than the 14±6 per meg observed at Kumukahi (20°N, 155°W), USA located at a similar latitude to HAT (Keeling et al., 2004).  surface temperature anomaly appeared. This correspondence is qualitatively reasonable since a decrease in the OHC change rate indicates a decrease in the net ocean heat uptake, and thus leading to an increase in surface temperature.
As mentioned above, the ratio of interannual variation to secular increase is larger for d(Ar/N2) than for OHC. To examine this difference, we estimated the interannual variation of the atmospheric d(Ar/N2) expected from the air-sea Ar and N2 fluxes caused by the interannual variation in the global OHC. We converted the change rate of OHC to that of d(Ar/N2) by assuming 190 a coefficient of 3.5x10 -23 or 3.0 x10 -23 per meg J -1 , which was derived from the following equations (Keeling et al., 1993).
Here, FAr (FN2) is the net sea-to-air Ar (N2) flux in moles m -2 s -1 , dCeq_Ar (dCeq_N2) is the temperature derivatives of the solubility of Ar (N2) in mole m -3 K -1 (Weiss, 1970), .̇ is the net air-to-sea heat flux in J m -2 , and CP is the heat capacity of seawater in J 195 m -3 K -1 . In eqs. (3) and (4), the assumption that surface waters fully equilibrate could lead to overestimation of the gas fluxes.
As a first approximation we modeled the global ocean as 1-box and assumed an average temperature of 3.5 or 7.5 °C. The temperatures of 3.5 and 7.5 °C, corresponding to the respective coefficients of 3.5x10 -23 (including deep water) and 3.0 x10 -23 per meg J -1 (not including deep water), were the average values of the ocean model shown in Fig. 1 of Bereiter et al. (2018) for the present-day ocean. We also assumed constant CP and salinity of 3.9x10 6 J m -3 K -1 and 35 ‰, respectively. To convert 200 FAr and FN2 to changes in atmospheric d(Ar/N2), we employed 5.124 x 10 21 g for the total atmospheric mass of dry air (Trenberth, 1981), 28.97 g mol -1 for the mean molecular weight of dry air, 3.6 x 10 14 m 2 for the surface ocean area of the earth, and respective fractions of 0.00933 and 0.7808 for Ar and N2 in the atmosphere.
As a result, the interannual variation in the change rates of the global OHC was estimated to drive only 10 % of that of the atmospheric d(Ar/N2). This discrepancy would be due to the fact that the troposphere, and not only the ocean, does not mix 205 perfectly on a timescale of a year, and that the surface-ocean temperature anomalies would be a large source of interannual variation on a yearly timescale for the observed d(Ar/N2) in the near-surface air. Therefore, in order to interpret the relationship between the interannual variations of OHC and d(Ar/N2) more quantitatively, we need to consider the local heat anomalies and atmospheric transport effects. In other words, a much longer-term d(Ar/N2) variation is less affected by the surface-ocean anomalies and atmospheric transport, but will include influence signature of the average temperature from the surface to the 210 deeper part of the ocean. In this regard, Suga et al. (2008) estimated that the renewal time of permanent pycnocline water in the North Pacific to be 2 -4 years for eastern subtropical mode water, 2 and 5 -9 years for the lighter and denser subtropical mode water, respectively, and 10 -20, 20 -30 and >60 years for the lighter, middle and denser central mode water, respectively. Therefore, it is expected that the secular trend of the d(Ar/N2) values obtained during 2012 -2019 (8-years) in the present study does not reflect the deep OHC change. 215 Figure 5 shows the d(Ar/N2) values observed at TKB and HAT, along with their best-fit curves and interannual variations from trends. Therefore, we applied another best-fit curve consisting of the fundamental and its first harmonics and a linear secular 220 trend to the observed data (black dotted lines), to extract the secular trends as a component of the linear trend of the best-fit curve (thick black solid lines). The secular trends were found to be 0.75±0.30 and 0.89±0.60 per meg yr -1 at TKB and HAT, respectively. The uncertainties for TKB (HAT) were evaluated, taking into account the respective uncertainties of ±0.11 and ±0.28 per meg yr -1 (±0.53 and ±0.29 per meg yr -1 ) for the regression and the long-term stability of the standard air. Comparable secular trends of 0.70±0.33 per meg yr -1 for TKB and 1.06±0.61 per meg yr -1 for HAT were also obtained from the difference 225 in the annual average value between 2012 and 2019 for each station. Although it is possible these changes observed over the study period might not represent a long-term trend, but is part of a large interannual variation. Nevertheless, it would of interest to see if it is possible to obtain a scientifically "meaningful" OHC change based on the secular d(Ar/N2) trend reported earlier.
Accordingly, we used the secular trend at TKB (0.75±0.30 per meg yr -1 ), because of the smaller uncertainty than that obtained at HAT, to propose a method to estimate the global OHC or BDC change based on the surface and the stratospheric secular 230 trends of d(Ar/N2). For this purpose, we simulate gravitational separation effect on the d(Ar/N2) in the whole atmosphere by using a 2-D model described below in subsection 3.2.

Simulation of gravitational separation and its effect on d(Ar/N2) at ground surface
As mentioned in the Introduction, we have reported observational results of gravitational separation in the stratosphere 235 (Ishidoya et al., 2008(Ishidoya et al., , 2013Sugawara et al., 2018). Those results showed that stratospheric d(Ar/N2) also decreases rapidly with increasing altitude above the tropopause. Also, not only are there large year-to-year variations in the difference between the stratospheric and tropospheric d(Ar/N2) values, year-to-year variations in stratospheric d(Ar/N2) are much larger than those observed in the troposphere (Ishidoya et al., 2013;. Therefore, we need to explore the possibility of tropospheric d(Ar/N2) variations caused by changes in the stratospheric gravitational separation that influence the whole 240 troposphere. We have compared the observed and simulated gravitational separation of atmospheric major components above the tropopause in previous studies (Ishidoya et al., 2013Sugawara et al., 2018;Belikov et al., 2019), but we did not consider gravitational separation in the troposphere. Therefore, in this study we updated the SOCRATES model to calculate variations in Ar and N2 from the surface to 120 km, taking into account molecular diffusion processes generating gravitational separation. Here, we describe only those modifications we made to the model for the calculation of d(Ar/N2). Additional 245 description of the SOCRATES model is presented in Appendix A.
First, our 2-D model was expanded to be able to calculate explicitly any inert gas components and their isotopic ratios and elemental ratios, including d(Ar/N2). In our previous model studies, we used 44 CO2 and 45 CO2 to reproduce gravitational separation for the sake of simplicity. However, it is now necessary to explicitly include the molecular diffusion coefficients of all gas components in order to reproduce the molecular diffusion processes in the stratosphere more accurately, because the 250 molecular diffusion coefficient is dependent on the molecular mass. In this study, the molecular diffusion coefficient of the gas component A in air (DA) was calculated using the following equation (Poling, et al., 2001). Second, a new δ value has been defined and used for d(Ar/N2) in this study. Usually, the isotopic ratios or elemental ratios of 260 the atmospheric major compositions are expressed as values relative to their ratios in the troposphere. Therefore, it is common to treat the d values at the ground surface as zero. However, high-precision analyses have recently revealed that, as is the case in d(Ar/N2), the d values at the ground surface are not always constant, and they have seasonal and inter-annual variations, as well as long-term trends (Keeling et al., 2004;Cassar et al., 2008;Blaine, 2005;Bent, 2014). The purpose of our analysis here is to evaluate the effect of gravitational separation on the d value at the ground surface using the 2-D model. Therefore, it is 265 no longer appropriate to assume that the d value at the ground surface is always zero. Observationally, small variations at the ground surface can be detected by using specific gas cylinders as constant references. In a numerical model, a constant reference is also needed for evaluating the effects of gravitational separation on the d value at the ground surface. In our 2-D model, we used the ratio of total amount (M) of each gas component in the model atmosphere as the constant reference and defined a new d value, δΩ, by the following equation. 270 The total amount was calculated by integrating from the ground surface to the altitude of 120 km. According to this definition, δΩ will be zero if there are neither molecular separations nor sinks/sources in the entire atmosphere. In the actual atmosphere, δΩ becomes negative in the stratosphere due to gravitational separation, but will be a small positive value in the troposphere.
A 50-year-long spin-up calculation was carried out for δ(Ar/N2)Ω under steady-state condition, and we found that δ(Ar/N2)Ω 275 in the troposphere reached a steady-state value of about 30 per meg. In other words, the tropospheric Ar/N2 ratio is enriched by about 30 per meg relative to the homogenous atmosphere due to the atmospheric gravitational separation. If the gravitational separation is strengthened due to atmospheric circulation changes in the stratosphere, δΩ at the ground surface will slightly increase. On the other hand, if the gravitational separation is weakened in the stratosphere, δΩ at the ground surface will decrease. Therefore, by using δΩ, it is possible to examine how a change in the gravitational separation in the stratosphere 280 affects δΩ at the ground surface. Meridional distribution of δ(Ar/N2)Ω calculated using the updated SOCRATES model, and its comparison with the stratospheric δ(Ar/N2) observed in our previous studies, are presented in Appendix B.
Third, age of air (AoA) is calculated using an idealised tracer. In our previous model calculations, the AoA was calculated from the CO2 mole fraction. However, since the actual increase in CO2 mole fraction given at the ground surface is not a linear increase but includes non-linear trends and inter-annual fluctuations, it is not possible to estimate the correct AoA from the 285 simple lag time method (e.g. Waugh and Hall, 2002). Therefore, in this study, we introduced an inert idealised tracer that increases linearly at the ground surface, and calculated the AoA in the stratosphere from the mole fraction of this idealised tracer.

Secular trends in the observed and simulated d(Ar/N2) and its implication for changes in OHC and BDC
To examine how d(Ar/N2)W is influenced by changes in the BDC, model simulations were made by arbitrarily changing the 290 mean meridional circulation (MMC) represented by mass stream function in the 2-D model so that AoA changed by ±0.02 yrs yr -1 at 35 km over the northern mid-latitudes. The simulated AoA values for the 35-km height are shown in Fig. 6. Here the positive and negative rates correspond to the weakening and enhancement of the BDC simulations, respectively. These simulations with fixed horizontal mixing are not enough to represent the mechanism to drive the observed AoA change in the real atmosphere, but it does serve as a kind of sensitivity test. In this regard, we present the annual mean meridional distribution 295 of the AoA trend for the weakened and enhanced BDC simulations in Appendix B. Based on some past studies, we have adopted a value of ±0.02 yrs yr -1 as AoA trends for the BDC simulations : for example, Diallo et al (2012) reported an increase in AoA of about 0.03 yrs yr -1 in the middle stratosphere for the period 1989-2010 based on the ERA-Interim reanalysis data, while Fritsch et al. (2016) reported AoA unchanged within uncertainties over the northern mid-latitudes for the period 1975-2016 based on the observations of CO2 and SF6 updated from Engel et al. (2009), andWaugh (2009) showed, based on results 300 from chemistry-climate models, a negative AoA trend between -0.005 and -0.02 yrs yr -1 for the time period and location considered by Engel et al. (2009). The secular AoA trends forced in our study fall within the range of these studies, although the observation periods do not overlap with each other. We regard these simulations to be a first step in investigating the effect of gravitational separation of the whole atmosphere on d(Ar/N2) at the surface. Secular changes in the simulated d(Ar/N2)W at the surface obtained from the weakened and enhanced BDC simulation, and 305 those at 35 km are also shown in Fig. 6. As can be seen from the figure, clear inverse relationships are found between the secular trends of d(Ar/N2)W at the surface and those at 35 km. In the weakened BDC simulation, d(Ar/N2)W changes secularly by 0.15 and -4.5 per meg yr -1 , respectively, at the surface and 35 km. In the enhanced BDC simulation, respective secular changes in d(Ar/N2)W by -0.13 and 4.0 per meg yr -1 at the surface and 35 km are found. As discussed above, atmospheric d(Ar/N2) is estimated to increase by 3.5x10 -23 per meg when 1 J of heat is added to a 3.5 °C ocean water mass. Based on this 310 relationship, when OHC increases by 93 ZJ during 2012 -2019, which is the global 0 -2000 m OHC increase from the NOAA/NCEI data, the increase rate in atmospheric d(Ar/N2) is estimated to be 0.46 per meg yr -1 . Therefore, 0.15 (-0.13) per meg yr -1 for the surface d(Ar/N2)W due to the weakening (enhancement) of the BDC is a non-negligible trend compared with the expected changes in atmospheric d(Ar/N2) due to the global OHC increase. Figure 7 shows the same annual average values and secular trend of atmospheric d(Ar/N2) at TKB as in Fig. 5. We calculated 315 an increase rate of d(Ar/N2)cor by subtracting 0.15 per meg yr -1 from the observed secular trend of 0.75±0.30 per meg yr -1 to remove the estimated effects of d(Ar/N2)W increase obtained from the weakened BDC simulation. Thus, the derived secular increase rate of d(Ar/N2)cor, which represents the increase rate of d(Ar/N2) driven only by the OHC change, is 0.60±0.30 per meg yr -1 . We can convert this increase rate to the global OHC increase rate by assuming above-mentioned coefficient of 3.5 (3.0) x10 -23 per meg J -1 (hereafter referred to as the increase rate of "OHCArN2"). The obtained secular increase of OHCArN2 by 320 17.1±8.6 (20.2±10) ZJ yr -1 is shown in Fig. 7, together with that of the global 0 -2000 m OHC reported by NOAA/NCEI (hereafter referred to as "OHCoc"). The increase rate of OHCArN2 is consistent with the OHCoc value of 12.2±1.2 ZJ yr -1 within uncertainties. Recently, by extracting a solubility-driven component of APO from their atmospheric O2/N2 ratio and CO2 measurements and ocean model simulations, Resplandy et al. (2019) estimated an increase rate of the global OHC to be 12.9±7.9 ZJ yr -1 for the period 1991 -2016 that is consistent, within the uncertainties, with the OHCoc during the same period. 325 The results of the present study and Resplandy et al. (2019) suggest the usefulness of atmospheric observations for independent confirmation of ocean heat uptake estimated from ocean temperature measurements.
As described above, the increase rate in OHCArN2 obtained by using d(Ar/N2)W from the weakening BDC simulation agrees with the OHCoc estimate. On the other hand, the increase rate in OHCArN2 of 25.1±8.6 (29.6±10) ZJ yr -1 , obtained by using d(Ar/N2)W from enhanced BDC simulation assuming the coefficient of 3.5 (3.0) x10 -23 per meg J -1 , is significantly larger than 330 the OHCoc estimate. In general, modeling studies have pointed to an accelerating BDC due to anthropogenic climate change (e.g. Austin and Li, 2006). However, the increase rate in OHCoc agrees well OHCArN2 based on the weakened BDC simulation rather than on the enhanced BDC simulation. In this regard, several balloon-borne observations (Ray et al., 2014) and the ERA-Interim reanalysis data (Diallo et al., 2012) have suggested aging of air in the northern hemispheric midlatitude midstratosphere, implying a slowdown in the deep northern hemispheric branch of the BDC. Garfinkel et al. (2017) analyzed a 335 series of chemistry-climate model experiments conducted for the period January 1960 through 2014. They found structural changes in BDC have occurred in the BDC since 1980s; BDC accelerated in the lower stratosphere in the northern hemisphere and tropics but not in the mid-stratosphere, and specifically since 1992, mean age increased by 0.12 year in the mid-stratosphere of the northern hemispheric mid-latitude and tropical mass upwelling has slowed down by 2 %. As discussed in our previous study , gravitational separation of the stratospheric air is sensitive to changes in mean meridional 340 circulation rather than horizontal mixing. Therefore, it is thought that the increase of mean age in the mid-stratosphere in the northern mid-latitude and the slowdown of tropical upwelling, suggested by the recent observational and modeling studies, are consistent with our weakened BDC simulation of the d(Ar/N2)W.
We also carried out an additional simulation of d(Ar/N2)W considering the interannual variation in the BDC, to examine its effect on the large interannual variations of the observed d(Ar/N2) seen in Fig. 4. As the BDC interannual variation, we simply 345 assume a 10 % change in the MMC intensity with a 3-year cycle by changing the mass stream function in the 2-D model. Flury et al. (2013) reported that the speeds of BDC towards the NH and SH show interannual variabilities with amplitudes of about 21 and 10 %, respectively, and that the amplitude of variability in the ascent rate at the Equator is 21 %. Therefore, the 10% change in the MMC in our simulation is not an unrealistic assumption. As a result, the increase rate of d(Ar/N2)W showed an interannual variation with a peak-to-peak amplitude of ±0.4 per meg yr -1 . This interannual variation is non-negligible but still 350 rather small compared to that found in the observed d(Ar/N2), about ±4.5 per meg yr -1 , as seen from Fig. 4. Therefore, it is suggested that the interannual variation of surface d(Ar/N2) is driven mainly by solubility change in seawater and/or that the assumed ±10 % changes in the MMC intensity in our model is smaller than the actual interannual variation of the MMC. In addition, analytical artifacts and the extremely difficult challenge of maintaining stable standard air for d(Ar/N2) could also be cause(s) of the interannual variation. Therefore, further studies are needed to achieve quantitative understanding of the 355 discrepancy.
Consequently, if we regard OHCoc as representing "true" global OHC, then we can estimate secular trend and interannual variations in the BDC from d(Ar/N2) observed at the surface. Conversely, if the simulated d(Ar/N2)W represents "true" BDC changes, then we can validate OHCoc by using OHCArN2. Therefore, the surface d(Ar/N2) is found to be a unique tracer for detecting changes in a spatiotemporally-integrated OHC and BDC. Not only precise observations of d(Ar/N2) over longer 360 periods of time, but also improvements in d(Ar/N2)W simulation using 3-D atmospheric transport models are important for future research advances in this field. For example, recently, Birner et al. (2020) simulated lower stratospheric d(Ar/N2) using the TOMCAT/SLIMCAT 3-D chemical transport model, which has been updated to include gravitational fractionation of gases. Such activities will contribute significantly to understanding the mechanisms of spatiotemporal variations in the stratospheric and the surface d(Ar/N2). 365

Conclusions
We have been carrying out systematic measurements of atmospheric d(Ar/N2) at TKB and HAT since 2012, COI and SYO since 2013 and 2016, respectively. Clear seasonal cycles of d(Ar/N2) with summertime maximum were found at TKB, COI and SYO, with peak-to-peak amplitudes of 21±10, 11±4, 5±10 and 32±9 per meg at COI, TKB, HAT and SYO, respectively, which are in general agreement with those observed at similar latitudinal sites by other investigators. The observed d(Ar/N2) 370 at TKB and HAT, after subtracting seasonal cycles and shorter-term variations, showed significant interannual variations and slight, but detectable, secularly increasing trends. The observed interannual variations correlated positively with that of the global OHC, suggesting a strong correlation between large-area air-sea heat flux and the long-term change in d(Ar/N2).
further studies are needed to interpret the discrepancy. 375 We improved the 2-D model we used in previous studies to calculate gravitational separation in order to evaluate its effect on the long-term change in the surface d(Ar/N2). We simulated weakened and enhanced BDC by arbitrarily changing the MMC (represented by mass stream function in the model), and obtained effects of 0.15 and -0.13 per meg yr -1 , respectively, on the secular trend of the surface d(Ar/N2). If we apply the correction for gravitational separation to the secular trend of d(Ar/N2) observed at TKB, then an average increase rate of 0.60±0.30 per meg yr -1 for d(Ar/N2) during 2012-2019 was obtained by 380 assuming a weakening BDC. By using a conversion factor of 3.5 (3.0) x10 -23 per meg J -1 by crudely assuming a 1-box ocean with a temperature of 3.5 (7.5) °C, we obtained an increase rate in the global OHC of 17.1±8.6 (20.2±10) ZJ yr -1 for the 8-year period, which is consistent with 12.2±1.2 ZJ yr -1 reported by NOAA/NCEI from the Argo float observations. On the other hand, the increase rate in the global OHC, obtained from the secular trend of d(Ar/N2) at TKB under enhanced BDC, was found to be significantly larger than that from the Argo float observations. These results indicate that the surface d(Ar/N2) is a 385 unique tracer for spatiotemporally-integrated OHC and BDC. Although an increase in global OHC is well known as an essential parameter to evaluate recent global warming, there is no method yet to adequately measure OHC via ocean temperature observations in the full-depth volume of the ocean. Long-term precise observations of d(Ar/N2) will meet this demand, after correction for gravitational separation of the whole atmosphere.

Appendix A: Additional description of the SOCRATES model
We performed numerical simulations using a 2-dimensional model of the middle atmosphere (SOCRATES) developed by the National Center for Atmospheric Research (NCAR) (Huang et al., 1998;Park et al., 1999;Khosravi et al., 2002). Since details of our model calculation for gravitational separation have already been described in Ishidoya et al. (2013) and Sugawara et al. (2018), only a brief description is given here. The model extends from the surface to 120 km altitude with a 1 km vertical 395 resolution, and from 85°S to 85°N latitude with a 5° latitudinal resolution. Time step of this model can be varied to accommodate chemical transport equations of various species. However, in our study, we used a time step of 1 day, because Ar and N2 molecules have no chemical reaction and the time constant of the molecular diffusion process is sufficiently longer than 1 day. The meridional and vertical eddy diffusivity coefficients (Kyy and Kzz) are parameterized by including three types of waves: planetary, gravity and tidal waves. The mass transport processes caused by molecular diffusion were originally taken 400 into account only above the mesosphere in SOCRATES, since the molecular diffusion effect was thought to be negligibly small in the troposphere and stratosphere, compared with the eddy diffusion effect. But since our study needed to include the effect of molecular diffusion in the lower atmosphere, however small compared to eddy diffusion, we simply lowered its vertical domain to the surface for the calculation of molecular diffusion. Molecular diffusion theory included in SOCRATES is based on Banks and Kockarts (1973).

Appendix B: Performance of the updated SOCRATES model
The annual mean meridional distribution of the d(Ar/N2)W value calculated using the updated SOCRATES model is shown in Fig. A1. The model shows that the vertical d(Ar/N2)W gradient increases with increasing height, especially at high latitude, which is consistent with the meridional distributions reported by Ishidoya et al. (2013) and Sugawara et al. (2018) calculated 410 using the original SOCRATES model. The basic structure of the vertical-meridional d(Ar/N2)W distribution can be interpreted as a result of a balance between the mass-dependent molecular diffusions and the mass-independent transport processes. Figure   A2 shows the vertical d(Ar/N2) profiles observed over Sanriku (39°N, 142°E), Japan on June 4, 2007 (Ishidoya et al., 2013) and Biak (1°S, 136°E), Indonesia on February 22-28, 2015 . The vertical d(Ar/N2)W profiles calculated by the updated SOCRATES model for the same seasons are also shown. As seen from the figure, the average vertical gradients 415 of the observed vertical d(Ar/N2) profiles are generally reproduced by the calculated d(Ar/N2)W, although the fluctuations in the observed profiles are not simulated. These fluctuations are also observed over Antarctica and could be attributed to the horizontal mixing of the stratospheric air over the region , their occurrences at low and middle latitudes are still not well understood. Overall however, we are satisfied with the performance of the updated SOCRATES model in reproducing the average vertical gradient of the d(Ar/N2) profiles from the surface to the stratosphere. 420 Figure A3 shows the annual mean meridional distribution of the AoA secular trend calculated using the updated SOCRATES model for weakened and enhanced BDC simulations. The simulations were made by arbitrarily changing the mass stream function in the model so that AoA changed by +0.02 (weakened BDC) or -0.02 (enhanced BDC) yrs yr -1 at 35 km over the northern mid-latitudes. The AoA trends are larger in the middle stratosphere than those in the lower stratosphere for the weakened BDC simulation. This altitudinal increase of the AoA trend is consistent with that derived from the ERA-Interim 425 (Diallo et al., 2012), however, the updated SOCRATES model cannot reproduce the negative AoA trend in the lower stratosphere found in Fig. 13 of Diallo et al. (2012). Ray et al. (2014) also reported increase and decrease of AoA in the middle and lower stratosphere, respectively. Therefore, simulations of the BDC secular trend using the updated SOCRATES model do not fully represent the details of the meridional distribution in the real atmosphere. Further studies are needed by using 3-D models for validating the d(Ar/N2)W secular trend discussed in the present study. 430 Data availability.
The observational data of d(Ar/N2) presented in this study can be accessed by contacting the corresponding author.

435
Author contributions.     In order to address the uncertainty of the secular trend, uncertainties around the regression and the long-term stability of the standard air (±1.6 per meg, blue circles in Fig. 1 (b)