What is the equivalent depth of the Pekeris mode?

Inspired by the detection of the Pekeris mode of atmospheric free oscillations by a recent 6 study, high-accuracy numerical calculations of the problem of determining the equivalent depth of atmospheric free oscillations are performed. Here, the computational method is 8 largely based on a previous study, but with modiﬁcations to improve the accuracy of the 9 calculation. Two equivalent depths are found, with values of 9.9 km and 6.6 km. The 10 former corresponds to the Lamb mode and the latter corresponds to the Pekeris mode. These values deviate from those obtained in the previous study, especially for the Pekeris 12 mode. The causes of this discrepancy is discussed, as well as the correspondence between the equivalent depths obtained in this study and that of the Pekeris mode detected in the 14 recent study.


Introduction
The explosive eruption of the Hunga Tonga-Hunga Ha'apai volcano on January 15, equivalent depth and the existence or non-existence of the Pekeris mode strongly depend 42 on the vertical temperature profile of the atmosphere, this discrepancy may be caused by 43 the difference between the temperature profile of USSA76 used by Salby and that of the 44 atmosphere on the day of the Tonga eruption when the Pekeris mode was detected by 45 Watanabe, et al. (2022). However, it is also possible that this discrepancy is due to an 46 accuracy problem in Salby (1979)'s calculation method itself, which is discussed in the next section. 48 In the present manuscript, we reexamine Salby (1979)'s calculation method of the 49 equivalent depth of the atmospheric free oscillations and propose a modifications to im-50 prove the accuracy of the calculation to determine the equivalent depth with higher ac-51 curacy. The remainder of the present paper is organized as follows. In Section 2, after HZ ′ − κZ = 0 (ζ = 0).
where z is the geometric altitude and H is a prescribed scale height. Note that H was 62 denoted byH in Salby (1979) but here we denote it by H to avoid confusion withH. 63 The parameter κ = (γ − 1)/γ, where γ is the specific heat ratio, and α = H/h, where h 64 is the equivalent depth. Also,H is a local scale height defined as, Here R 0 is the gas constant of the dry atmosphere, g 0 is the gravity acceleration, and T (ζ) 66 is the vertical temperature profile of the background field. These notations are changed 67 from Salby (1979) to be consistent with later descriptions in the present manuscript. 68 In Salby (1979), the vertical structure equation (1) and the boundary condition (2) 69 were not treated as they were, but the numerical calculation was done after applying the 70 following transform, and rewriting (1) and (2) as (equation (11) in Salby (1979)). Here, ξ is defined as, and k 2 (ζ; α) is a refractive index, which is defined as, Note that this explicit form was not written in Salby (1979). 75 titude, which seems to be ζ = 60 (not explicitly written in Salby (1979)), where the 77 radiation or evanescent boundary condition was imposed, and examined how well the 78 lower boundary condition (6) was satisfied with changing the value of α continuously.

79
There, the background temperature profile T (ζ) was set as described by USSA76, and 80 the prescribed scale height H was set as, where To overcome the above problems, we use only the following basic transform: Then, usingZ, the vertical structure equation (1) and the lower boundary condition (2) 113 can be written as, Now, by introducing (X, Y ) as we can derive the following equations: These are simultaneous ordinary differential equations for (X, Y ). Using (X, Y ), the lower 118 boundary condition (12) can be expressed as,
As the upper boundary condition, the radiation boundary condition or the evanescent 120 condition should be imposed. If we assume that it follows thatH =H t (constant) there. Then, from (14) and (15), we can derive the 122 following differential equation: Hence, introducing q as we can write the evanescent condition as if q < 0. Also, if q > 0, we can write the radiation condition as since we can choose one of the two solutions by choosing the sign of the frequency of 127 the disturbance without loss of generality. Because we are solving a linear homogeneous 128 problem, there is an arbitrariness of constant multiples in the solution. Therefore, we can set as, if q < 0. Also, if q > 0, we can set as, We can use either (23) or (24) as the starting condition at the point where ζ = ζ t , and 132 we can solve (16) and (17) in the decreasing direction of ζ. When (X, Y ) at ζ = 0 is 133 finally obtained, we can examine how well the lower boundary condition (18) is satisfied.

134
The calculation method introduced in this subsection will be referred to as the modified 135 Salby's method. there was still a problem that the gravity acceleration was assumed to be constant at 140 g 0 and the gas constant was treated as a constant (i.e., the mean molecular weight was 141 treated as a constant), even though the altitude range above 80 km was also treated.

142
In particular, the equivalent depth of the Pekeris mode may change if these effects are 143 taken into account. Hence, it is necessary to examine the case where these effects are 144 included. The vertical structure equation (1) and the lower boundary condition (2), 145 however, were derived with assuming that the gravitational acceleration and the gas 146 constant were constant in Salby (1979). Therefore, we will also perform the calculation 147 using the vertical structure equation and the lower boundary condition without these 148 assumptions. We begin with the following vertical structure equation and the lower 149 boundary condition in the log-pressure coordinate derived from the linearized primitive equations (Andrews, et al., 1987, equation (4.2.7a) and (4.2.7b)): Here,ẑ = −H ln(p(ζ)/p(0)) and p(ζ) is the vertical profile of background pressure. Note plitude of the disturbance in the log-pressure coordinate through the following equation: Here, R = R * /M and M is the mean molecular weight considering altitude dependence.

158
Note that the log-pressure buoyancy frequency differs from the usual buoyancy frequency.

159
In addition, note also that the definition of the squared log-pressure buoyancy frequency,

160
(27), is different from that in Andrews, et al. (1987). This form of definition is derived by 161 considering the altitude dependence of R. See equation (6.17.23) of Gill (1982) and the 162 explanation preceding it for details. Note, however, that the symbols are used differently.

163
Since the following relationship: holds betweenẑ and ζ (where, g(ζ) is the gravity acceleration considering altitude depen-165 dence), we can rewrite (28) as, if we introduceĤ as, Hence, (25) can be rewritten in ζ-coordinate as, Also, since we can rewrite N 2 * H 2 as, the vertical structure equation (31) can be rewritten as, The lower boundary condition (26) can be expressed in ζ-coordinate as, Furthermore, noting that R(0)T (0)/g 0 =Ĥ(0)H since g 0 = g(0), we can rewrite (34) as, Similarly as the previous subsection, if we introduce V as, we can rewrite (33) into the following simultaneous ordinary differential equations for Similarly as the previous subsection, the upper boundary condition can be imposed as 177 follows with assuming that T (ζ)/M (ζ) = (T /M ) t (constant) where ζ ≥ ζ t for (25) to be 178 a differential equation with constant coefficients. Then, if we introducer as we can impose the evanescent condition as, ifr < 0. Also, ifr > 0, the radiation condition can be imposed as, Considering that we obtain from (29) and (36), we can set as, ifr < 0. Also, ifr > 0, we can set as, Then, U H corresponds to Z except for constant multiples and we compare the profiles of 194 these variables in the next section.

196
3.1 Numerical results using the modified Salby's method 197 We integrate the simultaneous ordinary differential equations (16) and (17)  derivative of T that is necessary to computeH ′ is evaluated by the central difference as Here, we set ∆ζ = 10 m/H. This ∆ζ setting is also used for the ζ decrement in the 205

Runge-Kutta integration. Note that extrapolation based on the T definition in USSA76
206 is used to compute (47) at ζ = ζ t and ζ = 0. We have checked the dependence of the 207 following results on ∆ζ and sufficient convergence have been confirmed. 208 We repeat the integration with changing α continuously in the range of 0.5 ≤ α ≤ 1.5 209 and evaluate the left-hand side of (18). Figure  the same as those of the previous subsection but here we take the altitude dependence of 232 the gravity acceleration g(ζ) and the mean molecular weight M (ζ) described in USSA76 (Fig. 4).

259
Since the drawn profiles are for |W | in Fig. 6, which does not correspond to |Z| 260 shown in Fig. 3, we cannot compare these profiles directly. Instead of |W |, Fig. 7   From the results shown above, as long as the vertical temperature profile defined 284 USSA76 is used, we can conclude that the equivalent depths of Lamb mode and Pekeris 285 mode are 9.9 km and 6.6 km, respectively, with an accuracy of two significant digits if 286 the upper boundary for the computation is above the lower edge of the thermosphere. In 287 the present manuscript, the altitude dependence of the gravity acceleration and the mean 288 molecular weight are taken into account. In the thermosphere, the specific heat ratio γ 289 should also change with height, but this is not taken into account (In USSA76, only the 290 information that γ = 1.4 is given, and the altitude dependence of γ is not described).   . 9d), but the corresponding equivalent depth 327 for the spectral peak seemed to be larger than 6.1 km. This may imply that the long-term