Modulation Effects of Internal-Wave Evolution on Acoustic Modal Intensity Fluctuations in a Shallow-Water Waveguide

: Internal solitary waves evolving with time in shallow water are known to affect sound propagation signiﬁcantly. Unlike prior work studying the acoustic effects of individual internal-wave properties separately, this paper elucidates and evaluates the inﬂuence of a complete evolution process of internal waves on acoustic ﬁelds both theoretically and by the coupled ocean-acoustic simulation. Two evolving wave properties considered here are shape deformations including the variations of wave amplitudes and widths and packet dispersion manifested as the increasing wavelength (i.e., the distance between successive solitons). The acoustic modal intensity expressed by the Dyson series solution is reformulated to explicitly reveal the modulation effects induced by the deformation and dispersion of internal waves. Dispersion leads to modal interference and causes the intensity envelope to oscillate with the varying wavelength. Deformation modulates intensity in a non-oscillatory manner that is less predictable due to the complexity of amplitude and width variations. In the environment reconstructed from the ﬁeld observations of internal waves in the South China Sea, the modal intensity simulated by the parabolic-equation model exhibits pronounced modulation effects, where the modal interference due to dispersion dominates the intensity-envelope shape, and deformation affects the extremum positions of envelopes.


Introduction
Acoustic propagation in shallow water is susceptible to internal solitary waves (ISWs) [1] that can induce large perturbations into thermohaline structures and sound speeds [2][3][4]. From generation to dissipation, internal waves usually undergo a complex evolution process [5][6][7], in which shape deformations [8] and packet dispersion [9,10] are known to have substantial impacts on acoustic fields [11,12]. Here, the variations of wave shapes (e.g., amplitudes and characteristic widths) are termed deformations. Since individual solitons within the same packet typically propagate at different speeds [13], the wavelengths of internal waves change with time, which is called dispersion. Following the terminology used in physical oceanography, wavelength is defined as the distance between successive solitons [13].
Up to this point, the relationship between acoustic variability and individual properties of internal waves have been investigated extensively. The pioneer work conducted by Zhou et al. 1991 [11] attributed the observed anomalously large sound transmission loss to the resonant mode coupling that heavily depended on the wavelength of internal waves. Preisig and Duda (1997) [12] approximated a solitary wave by a square wave and explained the dependence of acoustic mode couplings on the width and amplitude of ISWs. Frank et al. 2004 [14] analyzed the sensitivity of the temporal intensity fluctuations of broadband airgun signals to the width and wavelength of ISWs. As a summary and extension of the previous work, Colosi (2008) [15] used the Dyson series solution [16] to the stochastic coupled mode equation [17,18] to derive a theoretical expression of acoustic modal intensity as a function of internal-wave parameters. The effects of amplitude, width, and wavelength on modal intensity were studied separately. Recently, Jiang et al. 2022 [19] considered the influence of irregular bathymetry on the propagation speed of ISWs in the environmental modeling and studied the fluctuations of sound intensity in the continental slope/shelf region. Gao et al. 2023 [20] investigated the acoustic temporal coherence in the presence of nonlinear internal waves and found that the correlation coefficients exhibited quasi-periodicity, the dominant period of which was the same as the periods of ISWs.
However, evolution actually involves the simultaneous changes of various properties of internal waves following physical rules, and the dependence of acoustic variability on the time-varying ISW parameters during a complete evolution process of internal waves remains to be studied. This paper features internal-wave evolution as an integral process rather than the isolated variations of individual wave properties. As a further refinement to the work in Colosi (2008) [15], we reformulate the Dyson-series representation of modal intensity to provide insights into the joint modulation effects of both shape deformation and packet dispersion on sound propagation. The different physical mechanisms by which deformation and dispersion affect the envelope of intensity fluctuations are analyzed, which are further confirmed in the modal intensity simulated by the parabolic-Equation (PE) model in an environment reproduced from an observed ISW packet in the South China Sea (SCS). In turn, the envelope variation of modal intensity obtained from acoustic data can be used to invert the dispersion rate of an ISW train.
The outline of this paper is as follows. In Section 2, we derive the expression of modal intensity as a function of internal-wave parameters using the Dyson series and elucidate the physical mechanisms by which deformation and dispersion modulate the intensity fluctuation. In Section 3, we describe the evolution process of an internal wave event observed in the SCS, and the sound-speed fields perturbed by ISWs are reconstructed for acoustic simulations. In Section 4, the sound fields during internal-wave evolution are calculated by the PE model, and the simulated modal intensity is analyzed to illustrate the modulation effect. In Section 5, we explain the modulation effects of deformation, dispersion, and both of them on modal intensity in the context of the Dyson-series representation obtained in Section 2. Finally, the conclusions of this paper are presented in Section 6.

Stochastic Coupled Mode Equation and Its Dyson Series Solution
Following the notations in Dozier and Tappert (1978a) [17] and Colosi (2008) [15], the Dyson-series representation of modal intensity in the single-scattering approximation is first derived in this section. The Helmholtz equation in the cylindrical coordinate system takes the form: where p(r, z) is the complex acoustic pressure,ρ(z) is the background density profile, ω is the acoustic angular frequency, and c(r, z) is the sound-speed field. Moreover, p(r, z) satisfies (1) the pressure-release boundary condition at the sea surface and (2) the continuity of pressure and vertical particle velocity at the water/bottom interface. Assuming that internal waves cause sound-speed perturbations δc(r, z) to the background profilec(z), the range-dependent c(r, z) is given by: The pressure p(r, z) in Equation (1) can be written in terms of normal modes, that is: where a n (r) and k n are the range-dependent modal coefficient and the real part of the horizontal wavenumber of the nth mode, respectively, N is the number of acoustic modes, and φ n (z) is the modal function in the background environment. Then, the rapidly oscillating phasor exp(il n r) is removed from a n (r) to get the demodulated modal coefficientã n (r), i.e., the complex envelope of a n (r):ã n (r) = a n (r) exp(−il n r), where l n = k n + iα n is the horizontal wavenumber of mode n and the imagery part α n accounting for sound attenuation can be estimated by the perturbation theory [21]. The range evolution ofã n (r) is governed by the stochastic coupled mode Equation [17,18]: where is the coupling matrix, c 0 is the reference sound speed, and k 0 = ω/c 0 is the reference wavenumber. Equation (5) does not have an exact analytical solution. We now approximateã n (r) by the Dyson series [16] and expandã n (r) at a receiver range r = r RX to the first order, yielding:ã whereã n (0) is the modal coefficient at the source range r = 0, and S (1) n,m is the first-order scattering matrix: where R ∧ n,m (k) in the scattering matrix is the Fourier transform of R n,m (r). The modal intensity |a n (r RX )| 2 is defined as the squared modulus of the modal coefficient a n (r RX ) [22]. Combining Equations (4) and (7) and keeping |a n (r RX )| 2 to the first order, the modal intensity in the Dyson-series form is written as: To obtain an analytical form of Equation (9), Colosi (2008) [15] used the Gaussian wave packet as the waveform of ISWs: where J is the total number of individual solitons within the packet and η(t; j), r p (t; j), and ∆(t; j) denote the wave amplitude, peak position, and characteristic width of soliton j varying with time t, respectively. The sound-speed perturbation δc(r, z, t) caused by internal waves is expressed as: where (dc/dz) ptnl is the potential sound-speed gradient [15,23] and W 1 (z) is the first-order internal-wave modal functionsatisfying the Taylor-Goldstein equation and the rigid-lid boundary condition [24]. Under the condition that the receiver range is large (r RX → ∞), and that the modal attenuation is small (α n , α m → 0), substituting Equation (11) into Equation (9) with some manipulations yields the analytical form of modal intensity, that is: where and

Modulation of Evolving ISWs on Modal Intensity Fluctuations
For clarity, the contribution from an arbitrary single incident mode m in Equation (12) is taken to discuss the influence of internal-wave evolution. The interference between the received mode n and the incident mode m is described by: and is further expanded in the ascending order as to the soliton numbering j: G n,m (t) = F n,m (t; 1) sin k m,n r p (t; 1) + F n,m (t; 2) sin k m,n r p (t; 2) + F n,m (t; 3) sin k m,n r p (t; 3) + · · · + F n,m (t; J) sin k m,n r p (t; J) .
Equation (16) provides information about the dependence of the mode-coupling strength on the wavelength of internal waves. That is, the fluctuation magnitude of G n,m (t) reaches to the maximum when the distance between successive solitons is equal to an integer multiple of the interference length L n,m determined by modes n and m, that is: which is called the resonant mode coupling [3,11,15]. However, the influence of packet dispersion on G n,m (t) remains implicit in Equation (16) and is not integrated with the effect of shape deformation reflected in F n,m (t; j). The modulation effects of deformation and dispersion on modal intensity can be made explicit by the following mathematical manipulations. Using trigonometric identities, the superposition of any two successive terms (e.g., j and j + 1) in Equation (16) can be written as: G n,m (t; j, j + 1) = F n,m (t; j) sin k m,n r p (t; j) + F n,m (t; j + 1) sin k m,n r p (t; j + 1) = Polychromatic term: poly n,m (t;j,j+1) 2F n,m (t; j) cos k m,n Λ p (t; j, j + 1) 2 sin[k m,n r c (t; j, j + 1)] where r c (t; j, j + 1) = r p (t; j) + r p (t; j + 1) 2 (19) is the position of the geometric center (also called the centroid) of the "dual-wave packet" composed of the two successive solitons j and j + 1, and is the distance between solitons, i.e., the wavelength of the dual-wave packet.
As formulated in Equation (18), G n,m (t; j, j + 1) is decomposed into the superposition of the polychromatic term: Carrier wave (21) and the monochromatic term: Carrier wave sin k m,n r p (t; j + 1) .
Equation (21) shows that poly n,m (t; j, j + 1) can be thought of as a high-frequency carrier wave with amplitude modulation. Typically, the centroid moves at speed dr c /dt of order 1 m/s, whereas the dispersion rate of a dual-wave train is dΛ p (t)/dt of much smaller order 0.1 m/s. Therefore, the movement of the centroid is responsible for the rapidly varying modal interference indicated by sin(k m,n r c (t)) in Equation (21).
Regarding the polychromatic envelope, there are two factors modulating the carrier wave together: one is dependent on packet dispersion; the other is associated with shape deformations. As the dual-wave packet disperses with time, the increasing Λ p (t) gives rise to the interference between modes n and m. Thus, the carrier wave is modulated in the form of harmonic oscillation indicated by cos(k m,n Λ p (t)/2) in Equation (21). In contrast to the dispersion typically manifested as the increase in wavelength, the variations of amplitude and width are heavily dependent on the environmental factors (e.g., stratification, topography, etc) of a specific sea area. Consequently, the modulation effect due to deformation may not be as regular as that caused by dispersion and is less predictable.
The composition of mono n,m (t; j, j + 1) in Equation (22) is simpler than that of poly n,m (t; j, j + 1). It is clear that only the movement of soliton j + 1 (i.e., the last wave in the dual-wave packet) leads to the modal interference at a single frequency, and that only deformation contributes to the amplitude modulation.
Applying the manipulation described by Equation (18) to Equation (16), we obtain another form of G n,m (t): G n,m (t) = poly n,m (t) + mono n,m (t), (23) where poly n,m (t) is the polychromatic component: and mono n,m (t) is the monochromatic component: By comparing Equation (16) with Equations (23)- (25), it is found that the use of trigonometric identities transforms G n,m (t) from the superposition form (see Equation (16)) into a form containing the product representation (see Equation (24)), thus accentuating the effect of dual-wave packets rather than individual solitons on acoustic intensity. Based on the above derivation, we next focus on the evolution process of a typical internal wave event observed in the SCS and investigate how shape deformations and packet dispersion modulate the acoustic intensity fluctuations in reality.

Experiment Description
The 2019 joint oceanographic-acoustic observation experiment was conducted on the continental shelf in the SCS. The mooring positions and the topography along the track are shown in Figure 1, and the bathymetry data used in this paper are from the GEBCO database [25]. There were a total of four moorings deployed along the track, and two of them were oceanographic moorings IW1 and IW2 equipped with temperature-pressure (TP) sensors, which were deployed 6.3 km apart. One was an acoustic mooring TX positioned 1.1 km away from IW1, on which a transducer was installed at 330.5-m depth and emitted five cycles of linear frequency-modulated (LFM) pulses sweeping from 400 to 500 Hz with a 10-s duration and a 20-s period every hour. The last one was a hybrid mooring RX equipped with both hydrophones and TP sensors and was deployed 18.9 km away from IW1.
The water depth decreased from 365 m (IW1) to 327 m (RX) with a mean value of 349 m and a mild slope of 0.1°. Internal waves propagated from RX to IW1, and sound transmission was in the opposite direction.

Oceanographic Observation
A total of 25 groups of internal waves were captured in the experiment, among which the fourth internal wave event (hereafter called Event 4) underwent a typical evolution process and is chosen for analysis in this paper. The temperature profiles measured at IW1, IW2, and RX during Event 4 are shown in Figure 2a-c, respectively, along with the 22.5°C isotherms denoted by blue lines. The individual solitons identified from data are annotated by red arrows. The in situ sound-speed profiles (SSPs) at IW1, IW2, and RX are displayed in Figure 2d-f, correspondingly.
The comparison of observations at IW1, IW2, and RX reveals that the internal waves exhibited shape deformations and packet dispersion as they propagated. The deformations were manifested as (1) the variations of wave amplitudes among moorings evidenced by the isothermal displacements and (2)

Reconstruction of Sound-Speed Fields
Before studying the sound propagation through internal-wave fields, we need to reconstruct the acoustic environment with high resolution in horizontal range and full coverage in vertical depth. Note that the time-independence of the modal coefficient a n (0) at range r = 0 in Equation (12)

Configuration and Method
The schematic of the waveguide model used in this paper is depicted in Figure 4 with simulation parameters and is in general agreement with the experimental environment. The continuous-wave (CW) signals at 450 Hz are emitted by a source near the seafloor of depth z b = 349 m and are received by a vertical line array (VLA) equipped with hydrophones at a distance of 17.8 km. Since a flat bottom is adopted in the simulation, the modeled source depth is scaled from the actual value z s,obs = 330.5 m to z s = z s,obs z b /z b,obs ≈ 320 m, where z b,obs = 360.5 m is the actual water depth at mooring TX. The bottom is modeled as a homogeneous halfspace with sound speed c b = 1750 m/s, density ρ b = 1.8 g/cm 3 , and attenuation rate α b = 0.5 dB/λ [27], where λ denotes the acoustic wavelength. With the sound speeds modeled by DEPTS being input parameters, acoustic fields p(r, z; t) are calculated by the computer code IFD [28] based on the PE model. Next, using the orthonormality relation of normal modes [23], mode decomposition is applied to the pressure p(r RX , z; t) received by the VLA at range r = r RX to obtain the modal intensity: where modal functions φ n (z) are calculated by the KRAKEN program [29].

Simulation Results and Analysis
A case study of the modulation effect in the first-order modal intensity |a 1 (r RX ; t)| 2 PE is carried out in this section. Figure 5a shows the |a 1 (r RX ; t)| 2 PE (blue) as a function of time (the bottom x-axis) and of the centroid position (the top x-axis). Here, the geometric center of the dual-wave packet composed of solitons #1 and #2 is designated as the centroid since these two solitons are the most prominent in the wave train, and the following solitons are not well-developed until 22:25:32 as shown in Figure 3d. The upper and lower envelopes are denoted by red-dashed lines. The peak-to-peak values (PPVs) defined as the difference between the upper and lower envelopes are used to quantify the fluctuation magnitude and are displayed by the green line on the right y-axis. It can be seen that |a 1 (r RX ; t)| 2 PE exhibits quasi-periodicity with a dominant period of about 30 min. The high-frequency oscillation is modulated by a slowly varying process as evidenced by the envelopes. The PPVs suggest that the fluctuation magnitude increases since internal waves enter the track at 19:11 and reaches its maximum at 20:25, after which it starts to decrease.
With regard to the range scale on the top x-axis in Figure 5a, it was found that after |a 1 (r RX ; t)| 2 PE oscillates at the dominant frequency for one cycle, the centroid moves about 2.65 km, close to the interference length L 1,2 = 2π/|k 1 − k 2 | = 2.84 km determined by the received mode n = 1 and the incident mode m = 2. The wavenumber spectrum of |a 1 (r RX ; t)| 2 PE is shown in Figure 5b with three pronounced peaks. The first spectral peak appears at k = 2.21 × 10 −2 rad/m, close to the horizontal wavenumber difference k 1,2 = k 1 − k 2 , and thus, represents the interference between modes n = 1 and m = 2, which is annotated by "Modes 1&2" in Figure 5b. Similarly, the other two peaks corresponding to the interference between modes 1&3 and between modes 1&4 are also identified, and it is clear that the first peak is at least 7 dB higher than others.
The above analyses in both the range and wavenumber domains demonstrate that the interference between modes 1 and 2 contributes most to the intensity fluctuation, and thus, the modulation effect found in the time domain. The formation mechanism of the highly variable fluctuation magnitude shown in Figure 5a will be investigated in the next section.

Study Approach
The following simplifications are made here to facilitate the theoretical interpretation of the modulation effect found in Figure 5. Since the interference between modes 1 and 2 dominates the intensity fluctuation, the Dyson-series representation in Equation (12) only including G 1,2 (t) is calculated to explain the physical mechanism behind the modulation effect in the above PE-simulation results. Besides, the entire ISW train observed in the experiment is approximated by a Gaussian wave packet consisting of the two most prominent leading waves (i.e., solitons #1 and #2 in Figure 3).
We now express the modal intensity |a 1 (r RX )| 2 in the partial-Dyson series form: where only the contribution from the incident mode 2 is included in the series. Referring to Equations (23)-(25), G 1,2 (t) in Equation (27) is rewritten as: where poly 1,2 (t) =  (29) is the polychromatic component, and is the monochromatic component. The wave parameters of Event 4 are output from the DEPTS model [26] and are used to calculate the Dyson-series representation given by Equation (27). The wave amplitudes η(t; j), characteristic widths ∆(t; j), and crest positions r p (t; j) of solitons j = 1, 2 are presented in Figure 6a-c, respectively. The centroid position r c (t; 1, 2) and the wavelength Λ p (t; 1, 2) are shown in Figure 6d. It is clear that as the wave packet propagates along the track, η(t; 2) keeps increasing while η(t; 1) starts to decrease after reaching its maximum at 21:36; ∆(t; 2) keeps decreasing while ∆(t; 1) are much more stable; soliton #1 moves faster than soliton #2, and hence, Λ p (t; 1, 2) keeps increasing. In addition, Figures 7a,b illustrate the potential sound-speed gradient (dc/dz) ptnl and the first-order modal function W 1 (z) of internal waves, respectively, used to calculate the sound-speed perturbation δc(r, z; t) given by Equation (11). Next, we first analyze the influence of deformation and dispersion separately and then discuss the joint modulation effect of these two factors on modal intensity from the perspective of a complete ISW evolution process observed in reality.

Modulation Effect Induced by Deformation
To study the influence of deformation exclusively, packet dispersion is neglected in this section by setting Λ p (t; 1, 2) to remain constant after 20:05 when soliton #2 has entered the track. The modal intensities calculated by the PE model (|a 1 (r RX ; t)| 2 PE ) and by the partial Dyson series (|a 1 (r RX ; t)| 2 Dyson:2 ) are shown in Figure 8a by the blue and red lines, respectively. The PPVs of |a 1 (r RX ; t)| 2 PE and |a 1 (r RX ; t)| 2 Dyson:2 are given in Figure 8b with the same color scheme. It can be seen from Figure 8a that the two solutions have similar dominant frequencies, whereas the temporal variations of envelopes are different. Figure 8b shows that the relative change in the PPV of the Dyson series solution after 20:05 is only 0.12/0.1 = 1.2, obviously smaller than that of the PE result (0.115/0.03 ≈ 3.8). Consequently, the fluctuation magnitude of |a 1 (r RX ; t)| 2 Dyson:2 is less variable than that of |a 1 (r RX ; t)| 2 PE , and the modulation effect caused only by deformation is less pronounced. To explain the above observations, the interactions among the individual terms annotated in Equations (29) and (30) are examined in more details. Figure 9a,c depict the envelope (purple-solid) and carrier wave (green) of poly 1,2 (t), respectively, where the deformation-dependent amplitude 2F 1,2 (t; 1) (purple-dashed) and the dispersiondependent interference cos(k 2,1 Λ p (t)/2) (orange) are also given. The envelope and carrier of mono 1,2 (t) are shown in Figure 9b,d, respectively. The superposition of poly 1,2 (t) (blue) and mono 1,2 (t) (yellow) and the resulting G 1,2 (t) (red) are illustrated in Figure 9e. The comparison between Figure 9a,b indicates that the magnitude of the monochromatic envelop is smaller than that of the polychromatic one. Therefore, the modulation effect in G 1,2 (t) is mainly determined by the polychromatic envelope, as shown in Figure 9e. It is clear from Figure 9a that since Λ p (t) no longer changes with time after 20:05, cos(k 2,1 Λ p (t)/2) stops oscillating and remains constant. Next, deformation starts to solely determine the polychromatic envelope that drops to its minimum at 21:40 and then increases with only a minor variation. Consequently, the modulation effect in |a 1 (r RX ; t)| 2 Dyson:2 governed by the polychromatic envelope is different from that in |a 1 (r RX ; t)| 2 PE . The above discussion reveals that the modulation effect induced by deformation only is not pronounced as compared to the PE result. Referring to Equation (13), the coupling between modes n and m is significant when the wave width satisfies the following condition [15]: However, in reality, the significance and manifestation of the modulation effect induced by shape deformations are less predictable. This is because remarkable discrepancies may occur in the variations of widths for different solitons in the packet. The condition given by Equation (31) is hard to hold for all solitons simultaneously. Additionally, as formulated in Equations (14) and (13), the variations of wave amplitudes η(t; j) during a long-term evolution process also introduce uncertainties into the modulation effect. Specific to our case here, Figure 10 shows the oscillation amplitudes F 1,2 (t; j) of the interference components G 1,2 (t; j) for solitons j = 1, 2. By comparing Figure 6a,b with Figure 10, it can be seen that η(t; 1) dominates F 1,2 (t; 1) more than ∆(t; 1) because ∆(t; 1) only has a small relative change in its value. In contrast, for soliton 2, both η(t; 2) and ∆(t; 2) change with time significantly so that both together affect F 1,2 (t; 2).  Figure 10. Oscillation amplitudes F 1,2 (t; j) of the interference components G 1,2 (t; j) for solitons j = 1 (blue) and 2 (red).

Modulation Effect Induced by Dispersion
To study the effect of dispersion exclusively, shape deformations are neglected in this section by setting η(t; j) and ∆(t; j) (j = 1, 2) to remain constant after 20:05. Figure 11a depicts the modal intensities calculated by the PE model and the Dyson series, and Figure 11b shows the PPVs of intensity fluctuations. As shown in Figure 11a, |a 1 (r RX ; t)| 2 PE and |a 1 (r RX ; t)| 2 Dyson:2 show good agreement in both the dominant frequency and the envelope variation. Figure 11b reveals the qualitative agreement in the trends of the fluctuation magnitudes of the two solutions. As compared to Figure 8b, the relative change in the PPVs of |a 1 (r RX ; t)| 2 Dyson:2 after 20:05 is 0.1/0.04 = 2.5 in Figure 11b, closer to that of the PE solution (0.115/0.03 ≈ 3.8), which leads to a more apparent modulation effect in |a 1 (r RX ; t)| 2 Dyson:2 , as shown in Figure 11a. The terms annotated in Equations (29) and (30) are presented in Figure 12 with the same format as Figure 9. The polychromatic envelope in Figure 12a is larger than the monochromatic one in Figure 12b, and thus, dominates the modulation effect in G 1,2 (t) and |a 1 (r RX ; t)| 2 Dyson:2 . η(t; j) and ∆(t; j) have stopped changing since 20:05, after which F 1,2 (t; 1), F 1,2 (t; 2), and F 1,2 (t; 2) − F 1,2 (t; 1) also remain constant, as indicated by the purpledashed line in Figure 12a and by the purple line in Figure 12b. Then, the modal interference caused by packet dispersion starts to solely control the polychromatic envelope. The increment in Λ p (t) between 20:05 and 22:25 is less than half of the interference length L 1,2 , and cos(k 2,1 Λ p (t)) increases monotonically during this period with a notable relative change larger than (−1)/(−0.4) = 2.5. The polychromatic envelope is similar to the temporal trend of the fluctuation magnitude of |a 1 (r RX ; t)| 2 PE given in Figure 11b, which reveals the dominant role of the dispersion-dependent interference in modulating the intensity fluctuation. (e) poly 1,2 (t) (blue), mono 1,2 (t) (yellow), and G 1,2 (t; 1, 2) (red).

Modulation Effects Induced by Both Deformation and Dispersion
In this section, both shape deformations and packet dispersion shown in Figure 6 are considered when calculating the Dyson-series representation of modal intensity denoted by Equation (27). The modal intensities calculated by the PE model and by the Dyson series are shown in Figure 13a, and the PPVs are given in Figure 13b, both of which are in the same format as Figures 8 and 11. Figure 13a shows the qualitative agreement in both the oscillation patterns and the envelope variations between the numerical and theoretical results. Figure 13b suggests that the relative change in the PPVs of |a 1 (r RX ; t)| 2 Dyson:2 after 20:05 is 0.105/0.04 = 2.63, close to that of |a 1 (r RX ; t)| 2 PE (0.115/0.03 ≈ 3.8). The individual terms annotated in Equations (29) and (30) are shown in Figure 14. Furthermore, its enlargement is shown with the orange and purple arrows in Figure 14a, denoting the trough positions of the dispersion-dependent interference and the polychromatic envelope, respectively. Similar to Figures 9 and 12, the polychromatic envelope here still dominates the modulation effects in G 1,2 (t) and |a 1 (r RX ; t)| 2 Dyson:2 . However, unlike the above two cases, the zoom plot in Figure 14a reveals that the minimum of the polychromatic envelope is delayed backward by about 30 min compared to that of cos(k 2,1 Λ p (t)). This is because 2F 1,2 (t; 1) contributes to the polychromatic envelope at the same time by multiplying cos(k 2,1 Λ p (t)) and then shifts its trough position originally at 20:00. It can be seen that the modal interference resulting from dispersion governs the envelope shape of the intensity fluctuation, while deformation affects some fine structures in the envelope, such as extremum positions.   (e) poly 1,2 (t) (blue), mono 1,2 (t) (yellow), and G 1,2 (t; 1, 2) (red). Shown also in Figure 14a is an enlargement, with the orange and purple arrows denoting the trough positions of the dispersion-dependent interference and the polychromatic envelope, respectively.

Conclusions and Discussions
An analytical expression of modal intensity is reformulated based on the Dyson-series representation in this paper to describe the joint modulation effects of the shape deformation and packet dispersion of internal waves on sound propagation. Dispersion modulates modal intensity fluctuations by causing the interference between acoustic modes so that the intensity envelope exhibits a regular pattern of harmonic oscillation. Deformation modulates sound intensity in a non-oscillating manner, and the envelope may be irregular and less predictable due to the difficulty in fully understanding the realistic evolution process of internal-wave amplitude and width.
A coupled ocean-acoustic simulation is performed by first reconstructing the soundspeed fields from the experimental observations in the SCS and then modeling sound propagation using the PE model. The simulated modal intensity oscillates in a quasi-periodical pattern, and the dominant frequency is related to the interference length determined by the two modes with the most pronounced coupling strength. The modulation effect is manifested as the highly variable envelope of intensity fluctuations. The dispersion of ISW packets causes the periodical modal interference and governs the envelope shape, whereas internal-wave deformation is less irregular and affects the extremum positions of the envelope.
Due to the insufficient transmission of acoustic signals in the experiment (i.e., only one emission containing five pulses per hour), both the modal intensity fluctuation and the modulation effect are hard to detect in the acoustic data. A direct model/data comparison is currently unavailable and is not the objective of this study. Nevertheless, the high-fidelity sound-speed fields reconstructed from the experimental observations make the acoustic simulation not impractical. This paper provides a priori physical information for the target detection in dynamic marine environments and shows the possibility of inverting the internal-wave dispersion process from the envelope of acoustic intensity fluctuations. A joint oceanographic-acoustic experiment with intensive transmissions of sound pulses is still needed in the future to investigate the fine structure of the fluctuating acoustic field caused by internal waves.

Acknowledgments:
The authors would like to thank all participants working for the experiment. Their great efforts contributed to the valuable data for this paper. The authors also thank the GEBCO Compilation Group (2023) GEBCO 2023 Grid (accessed on 15 April 2023) (doi:10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b) for providing the bathymetry data.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: