[Note: In "Author's response", figure and equation numbers are those in the original manuscript to make them understandable without the marked revised manuscript. In "Author's changes in manuscript", figure, equation, and line numbers generally correspond to those in the marked revised manuscript, so that the referees and the editor can check the corresponding changes.]

-------------------------

<Referee comments>

As part II of a series of manuscripts, this one applies the linear stochastic model proposed in Part II to ocean data to identify sources of non-harmonic internal tides. It is a very good idea to use stochastic modelling to analyze measured data, and since this model is linear, the inversion problem is possible. 

My major concern is the validity of this linear modelling applied to a highly nonlinear system. The author did not justify the model in a controlled experiment. Maybe this is done in Part I? e.g., the impact of the background eddies, which vary in space and time, modify the Matrix A in (25). Whether a simple linear first-order SDE can capture this effect?

<Author's response>

Thank you for your comments, and I apologize that I did not justify the linear approach. I added a new appendix to compare the order of magnitude of the phase-speed modulation caused by mesoscale variability and the nonlinear triad wave interaction. The result suggests that the nonlinear triad wave interaction is an order of magnitude of smaller than the modulation effect for VM1 semidiurnal internal tides, both in the deep ocean (within the modelled region) and the PIL200 location on the continental shelf. This provides justification for using a combination of linear models as a first approximation.

I would also point out that your argument is somewhat misleading, although I agree that the mesoscale variability is highly nonlinear. This is because the cumulative effects of wave modulation caused by strongly nonlinear processes are not necessarily nonlinear. This is known in the study field called "wave propagation in random media". For example, turbulence and short stochastic internal waves (approximately represented by the well-known Garrett-Munk spectrum) are nonlinear, but an acoustic signal modulated by these processes can be modelled well by linear methods (see e.g., Colosi 2016). Of course, internal tides are not as linear as acoustic waves, but even large-amplitude internal solitary waves and high-frequency internal waves near the buoyancy frequency are often weakly nonlinear, except on shallow shelves and upon wave breaking.

In my SDE formulation, the matrix A represents the background state without eddies, and the effects of eddies are aggregated in c' (eddy-modulation of the phase speed) in Eq. (27). When the observed internal tides travel over many eddies and/or the observed signal consists of internal tides from many independent sources, statistical principles (the central limit theorem) make the observation insensitive to the details of phase deviation and c' other than their variance. Furthermore, as a by-product of the analysis in the new appendix, the derivation yielded Eq. (27), which provides another justification of the linear first-order SDE approach. So, I believe the simple linear first-order SDE is a reasonable first approximation when the nonlinear triad wave interaction is negligible compared to the phase-speed modulation effect. 

<Author's changes in manuscript>

- New Appendix A was added to compare the order of magnitude of the phase-speed modulation and nonlinear triad wave interaction.

- In new Appendix A, it is mentioned that the cumulative effects of wave modulation caused by strongly nonlinear processes are not necessarily nonlinear in general. 

- l.1225: "Note that the matrices A and B are calculated for the background conditions, and that c' aggregates the effects of interannual and mesoscale variabilities. The processes inducing c' can be strongly nonlinear, but the wave modulation process under given c' is approximately linear, as shown in Appendix A."

-------------------------

<Referee comments>

Section 3: A suggestion is to reduce this section to provide only the necessary equations that are used in analyzing data. 

<Author's response>

Thank you for your suggestion. Considering your and another referee's comments, I moved the detailed points in Section 3.1 and most of the derivation in Sections 3.2-3.5 to new Appendices C, D, and E.

<Author's changes in manuscript>

- The details regarding the use of R^{1/2} in Section 3.1 was moved to new Appendix C.

- Most of the derivation in Section 3.2 was moved to new Appendix D.

- The original Section 3.3 and the derivation in Sections 3.4-3.5 were moved to new Appendix E.

- Corresponding changes were made to the text.

-------------------------

<Referee comments>

Also, many assumptions leading to the model are missing. e.g., in the definition of (3), \Theta' is not a small perturbation of \phi. Does this impact the modelling?

<Author's response>

No, large \Theta' does not impact the modelling. This is essential because the variance of \Theta' (phase spread) keeps growing with wave propagation. The wrapped normal distribution, Eq. (2), is used in this study because it can handle arbitrary large phase spread.

Regarding missing assumptions, I have added new Appendix A to show that nonlinear effects are negligible as a first approximation, and the derivation provides assumptions leading to the use of the linear models. Overall, the important assumptions are as follows.

- Mesoscale and inter-annual variability, which induce phase-speed variability, are slowly varying in time compared to the wave variability.

- The variability of celerity is small compared to the background value. (This was assumed in the stochastic phase modelling in the original manuscript.)

- Nonlinear triad wave interaction is negligible compared to the wave modulation due to phase-speed variability.

<Author's changes in manuscript>

- New Appendix A was added to show assumptions leading to the use of the linear models.

- l.359: "Note that P_{\theta\theta} can grow without a limit, but this does not cause any problem because the wrapped normal distribution, Eq. (2), can be used with arbitrary large phase spread \sigma_j."

-------------------------

<Referee comments>

(15) is strange; to reach a statistically steady state, damping in this system is required. Is there a damping effect in A?

<Author's response>

I am confused by your comment, because there is no A in Eq. (15).

If your comment is about Eq. (15), damping is included in L (see Appendix A in the original manuscript). Even without damping, a forced oscillatory system can reach a stationary state with a constant amplitude unless forcing frequency is close to one of the natural frequencies of the system (resonance).

If your comment is about Eq. (25), it does not have to reach a stationary state (and many solutions of stochastic differential equations do not). For example, in the well-known random walk, the person can move further away from the initial position with time (or the total number of steps). In the case of internal tides, they become more random as they propagate through spatially and temporally varying ocean. The term "non-stationary" internal tides reflect this fact.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

However, dissipation in the primitive equation happens at small scales, while the forcing appears at large scales, then there must be a nonlinear effect to transfer energy across scales, which is missing in the current linear model.

<Author's response>

Thank you for your comment. Although the nonlinear energy cascade from internal tides to mixing is an important topic in oceanography, it is not important for this study as shown in new Appendix A. The hydrodynamic model in this paper includes only bottom friction, which is important for determining wave amplitudes on continental shelves. For bottom friction, the energy cascade occurs in the bottom boundary layer, which is parameterized in the vertical-mode formulation used in this paper. Then, the momentum (energy) loss is transferred to internal tides within the ocean interior through the Ekman transport and the unsteady version of the well-known spin-down of quasi-geostrophic flows through "inviscid" processes (e.g., Shimizu and Imberger 2009). So, the nonlinear cascade does not have to be modelled explicitly.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

(18a) and (21) look inconsistent.

<Author's response>

Thank you for your comment. It may not be intuitive, and that is why I show Eq. (21). If you formally apply the Fourier integral to (18a) for t = t_j to -\infty (because the adjoint model runs backwards in time), and use integration by parts to the LHS and the initial condition (18b), you will see that you get an equation of the form of Eq. (21). However, this involves some detailed points, and it is much easier to assume periodic motion in the original equation, Eq. (15), and derive the adjoint model from there.

<Author's changes in manuscript>

l.1147: "This may appear inconsistent with Eq. (D4), but can also be obtained by considering the Fourier integral of Eq. (D4a), and applying integration by parts to the left-hand-side and the "initial" condition Eq. (D4b), assuming \lambda = 0 for t > t_j."

-------------------------

<Referee comments>

Line 332: the relation d \theta = ..., should dispersion enter this relation?

<Author's response>

I am not sure if I understand your comment. If you refer to the wave dispersion of internal tides (e.g., caused by the Coriolis effects), it does not appear in the relation, because it is included in the baseline (unperturbed) solution (shown in new Appendix A). The relation you refer to is a simplified assumption made by Zaron and Egbert (2014), which has been supported by previous studies and further supported by the analysis in new Appendix A.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

Line 342, (27) is not an equation.

<Author's response>

I do not understand your intention. It may not appear as a standard differential equation, but it is a valid stochastic differential equation (see textbooks of stochastic differential equations, such as Sarkka and Solin 2019). Note that \theta'' and c' are stochastic variables.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

(32): P_\theta\theta tends to infinity as t tends to infinity, which looks unreasonable.

<Author's response>

Thank you for your comment. P_\theta\theta does not have to remain finite for very large t. For example, in the well-known random walk, the person can move further and further away from the initial position with time (or the total number of steps), so the variance of the person's position keeps increasing with time without a limit. In the case of internal tides, they become more random as they propagate through spatially and temporally varying ocean, if internal tides are not dampened by other processes. Very large P_\theta\theta does not cause a problem in the final result, because the variance of each wave component under the wrapped normal distribution (Eq. (4b)) reaches an upper limit as P_\theta\theta =\sigma_j^2 increases.

<Author's changes in manuscript>

- l.359: "Note that P_{\theta\theta} can grow without a limit, but this does not cause any problem because the wrapped normal distribution, Eq. (2), can be used with arbitrary large phase spread \sigma_j."

-------------------------

Reference

Shimizu, K. and Imberger, J.: Damping mechanisms of internal waves in continuously stratified rotating basins, J. Fluid Mech., 637, 137–172, doi:10.1017/S0022112009008039, 2009.

[Note: In "Author's response", figure and equation numbers are those in the original manuscript to make them understandable without the marked revised manuscript. In "Author's changes in manuscript", figure, equation, and line numbers generally correspond to those in the marked revised manuscript, so that the referees and the editor can check the corresponding changes.]

Review of "Process-based modelling of nonharmonic internal tides using adjoint, statistical, and stochastic approaches. Part II: adjoint frequency response analysis, stochastic models, and synthesis" by Shimizu

---------- General Comments ----------

<Referee comments>

The author has developed an original approach to analysis of tidal variability and applied it to understand observations from the Australian shelf. He uses the model to quantify sources---source regions, source strengths, vertical mode, and frequency---of tidal variability observed at a mooring. Although the results are very specific to this site and the regional setting, the overall approach is, in principle, more generally applicable. Furthermore, the approach has enabled a reasonably thorough qualitative discussion of the mechanisms of refraction along the internal propagation paths, which includes a nice discussion and analysis of the sensitivity to many of the simplifying assumptions. Although I suppose there is a narrow audience of readers with interest in this type of detailed analysis, the creativity and novelty of the approach should inspire substantial follow-on work. For this reason I think it should be published.

As I understand it, the basic components of described approach are as follows: (1) There is a deterministic model for the internal tides, which assumes linear dynamics (based on the author's previous work, and outlined in the Appendix). This model, and its adjoint, are used to compute the sensitivity of the waves at the observation site to the distributed sources surrounding it throughout the nearby Eastern Indian Ocean. The (mean) wave propagation properties are assumed to be constant in time. (2) There is a model for the second-order statistics of the non-harmonic tide phase and phase-speed modulations, and their spatial correlations. This model, and its adjoint, are integrated along ray-paths to describe (map) the statistics of the non-harmonic phase of the tide as it propagates from source regions (backward along the rays). (3) There is a model (developed in detail in Part I) which is used to relate the phase statistics of the individual wave sources to the statistics of the sum of the waves (at the observation site).

This manuscript (Part II) derives in detail components (1) and (2) above. Then it estimates the necessary input parameters from first principles (the barotropic-to-baroclinic tide forcing) and observations (phase speed variance and correlation scales), and then it proceeds to compare with the observed tidal variability and its sensitivity to the modeling assumptions (principally, the spatial correlation structure of the phase speed variability).

Overall, I found the manuscript too long, and rather hard to follow. Due to the complexity and multi-step development, it is essential for the author to reduce the verbage to the minimum necessary to communicate clearly. As it is, there are several sets of comments and caveats mentioned, which, while appropriately nuanced and apparently relevant, make it hard to follow the thread of the essential analysis.

It seems to me that the use of the adjoint model for the wave linear dyanmics (used to produce Fig 5) is relatively well-worn in the oceanography literature. The modeling of the along- and across-path covariance structure of the phase modulations (eqn's (30)-(38)) seems to be totally new at this level of precision and detail. I would suggest that the author consider breaking this Part II manuscript into two smaller pieces, one focused on the ray-tracing and modeling of the phase covariance, and then the other on using this covariance with the adjoint wave model to explain the observations. I think this manuscript is skillfully using a lot of innovative ideas, and it will have more impact if the it is broken down into simpler and more digestible pieces. Of course, this type of re-organization of the presentation will take considerable work.

Given the relatively narrow audience, maybe it is not worth the effort, but I hope the author will consider it.

In any case, I suggest this article be accepted after significant revisions to address my detailed comments, below. Some of my questions are answered by text, later in the manuscript, and reflect my misunderstandings. Nonetheless, I hope my comments will inform the author of the reaction of an interested reader, and guide him in making the manuscript more comprehensible.

<Author's response>

Thank you very much for your thorough reading and many constructive comments and suggestions. Also, thank you for recognizing the creativity and novelty of the work despite difficulty in reading.

Following your and another referee's comments, I made substantial changes to make the manuscript easier to read. This included moving most of the derivation and detailed points in Section 3.1-3.5 to new appendices, and removing non-essential details as much as possible. A short list of the changes are given below. As a result, the main body of the paper is about 20% shorter than the original manuscript.

I did consider splitting the whole study into three parts before submission. However, I had a difficulty that one of them does not have sufficient results as a stand-alone journal article in oceanography. For example, if I split the phase modelling part as your suggestion, none of the conclusions could be transferred to that paper. So, it would require time and effort to do new independent modelling and/or data analysis, and to write essentially another journal article. Your suggestion does make sense, and I would be interested in doing so if I had time and budget for it. However, for practical reasons, I decided not to split Part II further.

<Author's changes in manuscript>

- The length of Section 2 was halved.

- The details regarding the use of R^{1/2} in Section 3.1 was moved to new Appendix C.

- Most of the derivation in Section 3.2 was moved to new Appendix D.

- The original Section 3.3 and the derivation in Sections 3.4-3.5 were moved to new Appendix E.

- All the materials related to PDF and degrees of freedom were deleted. This includes the whole Section 4.6, the last paragraph in Section 5.6, and panels (c) and (d) in Fig. 10 in the original manuscript (now Fig. 11).

- The whole manuscript was edited thoroughly to remove non-essential details as much as possible, and to improve readability. As a result, Fig. 4 in the original manuscript was deleted.

---------- Detailed Comments ----------

<Referee comments>

The abstract says that a map of non-harmonic internal waves sources is identified from data on the Australian North West Shelf. The abstract could be clearer about exactly what kind of data are used.

<Author's response>

Thank you for your comment. I added the following sentence in Abstract.

<Author's changes in manuscript>

l.11: "Essential inputs of the model suite are barotropic tidal currents, background stratification, and the variance and spatial correlation of internal-tide phase speed."

-------------------------

<Referee comments>

l42: Does this sentence make sense? It is comparing "generation" with "amplitude modulation" and "phase modulation". The "generation" is of a different category than modulation.

<Author's response>

Thank you for your comment. I have revised the sentence as follows.

<Author's changes in manuscript>

l.44: "Although the variability of internal-tide generation can be substantial (Kerry et al., 2016), the amplitude variability is overall considered to be less important than the phase variability (Colosi and Munk, 2006, Zaron and Egbert, 2014)."

-------------------------

<Referee comments>

l97-l123: This is a long overview, but I don't feel like it has provided me with specifics needed to understand what is to come. Maybe it can be shortened or omitted.

<Author's response>

Thank you for your suggestion. I have halved the length of Section 2.

<Author's changes in manuscript>

- Section 2: "An overview of the proposed modelling framework is shown in Fig. 1. The key component is the statistical model developed in Part I. It calculates the statistics of nonharmonic internal tides by randomizing the phases (and optionally amplitudes) of individual internal-tide components arriving at an observation location from deterministic sources. For realistic oceanic applications, horizontal distributions of the sources and phase statistics are necessary. The source distribution can be modelled using an adjoint sensitivity model and barotropic tidal forcing. The implementation in this study uses a combination of numerical adjoint sensitivity modelling and the frequency response analysis from Fourier theory, referred to as "adjoint frequency response analysis". Currently, there appears to be no standard method to model the distribution of phase statistics. Since phase statistics vary with wave propagation (i.e., nonstationary), its process-based modelling appears to require a stochastic approach. The implementation in this study uses two stochastic models to model the spread of wave phases and the horizontal (two-dimensional) correlation of phase modulation, both of which are assumed to be caused by random variability of the phase speed. The final result is the statistics of nonharmonic internal tides, such as their PDFs (not shown in this paper) and the horizontally distributed sources of their variance."

- Caption of Fig. 1: "Overview of proposed modelling framework and its implementation in this study. Entire process applies two "filters": (i) to transform global and deterministic forcing from barotropic to individual baroclinic modes (forcing function), to forcing relevant only to a particular observation location (source function); and then (ii) to transform this forcing to response relevant only to random component of internal tides (nonharmonic variance source function)."

- Fig. 1 was updated.

-------------------------

<Referee comments>

Eqn (4a)-(4d): Where are these properties of the wrapped normal proven?

<Author's response>

Thank you for your comment. The reference to Part I was added.

<Author's changes in manuscript>

l.162: "Assuming tentatively that \sigma_j in Eq. (2) are known, and that all the wave components are independent, the expectation and variance of the complex-valued random amplitudes a_j e^{-i(\varphi_j+\Theta_j)} are (see Part I) ..."

-------------------------

<Referee comments>

Eqn (4b) and Eqn (5): This is confusing. At line 137, it says that A_j is deterministic. Doesn't this mean that A'_j is identically zero?

<Author's response>

No, A'_j is random because the phase \Theta_j is random, and because A_j' is the magnitude of the deviation from the mean of a_j e^{-i\Theta_j} on the complex plane. In response to this and other comments below, I added a new figure showing the wrapped normal distribution and the relationship among (a_j,\Theta_j), (r,\phi_j), and (A_j',\Theta_j'). I hope this makes clear why phase randomness makes A_j' a random variable, even when a_j is deterministic. Also, to clarify that A_j = a_j throughout the manuscript, A_j was replaced by a_j except Eq. (1).

<Author's changes in manuscript>

- New Figure 2 was added.

- A_j was replaced by a_j throughout the manuscript except Eq. (1).

- l.160: "Note also that \Theta' and \Theta_j' are random variables with zero mean unlike Part I, and that A, A', and A_j' are random variables even though a_j is deterministic (see Fig. 2 and Part I)."

-------------------------

<Referee comments>

l171: I am not understanding this. I don't understand what distinguishes \Theta'_j and \Theta''_j.

<Author's response>

I am sorry that this was unclear. In response to this and other comments, I added a new figure showing the wrapped normal distribution and the relationship among (a_j,\Theta_j), (r,\phi_j), and (A_j',\Theta_j'). To make the notation simpler, I also changed the definition of \Theta, so that the new \Theta is the same as the old \Theta''.

In addition, the explanation about Eq. (6) was revised, because it was not accurate in the original manuscript.

<Author's changes in manuscript>

- New Figure 2 was added.

- The explanation about Eq. (6) was revised.

- The mean phases were subtracted from \Theta, \Theta_j, \Theta', and \Theta_j', and \Theta_j'' was replaced by \Theta_j, throughout the manuscript.

- l.147: "Unlike Part I, the mean phase lags are subtracted from the total phase lags to make \Theta and \Theta_j random variables with zero mean, and only deterministic amplitudes A_j = a_j are hereafter considered for individual wave components."

- l.160: "Note also that \Theta' and \Theta_j' are random variables with zero mean unlike Part I, and that A, A', and A_j' are random variables even though a_j is deterministic (see Fig. 2 and Part I)."

------------------------

<Referee comments>

l195: n is a complex random number, right? I am not knowledgeable enough about the properties of complex random variables to be certain that the Cholesky-like decomposition mentioned at this line exists. Is R^{1/2} always defined for complex n? Is it complex-valued or real-valued?

<Author's response>

Yes, n is a complex-valued vector. However, R is real-valued, whose (i,j) components are given by Eq. (7). It is not intuitive, but Eq. (8) shows that the expectation of the product of n and its complex conjugate is real-valued. To derive this convenient relationship, it is essential that \Delta\Theta'' is defined with zero mean.

<Author's changes in manuscript>

- l.192: "Note that this relationship makes the correlation coefficients R_{ij} real-valued, although the original variable, e^{-i\Delta\Theta}, is complex-valued. To derive this convenient relationship, the definition of \Theta_j is changed from Part I to have zero mean." 

- l.212: "Note that n is complex-valued, but R is real-valued because of Eq. (8)."

-------------------------

<Referee comments>

l206-l220: This section points out the non-uniqueness of the matrix square root. Is R always real-valued?

<Author's response>

Yes, R is always real-valued, as explained in my response to your last comment.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

l221-l227: I'm afraid I don't understand this..

<Author's response>

I am sorry that the section was unclear. I revised the paragraph as follows. I hope it is understandable now. (The detailed points regarding the treatment of horizontal correlation using R^{1/2} were moved to new Appendix C.)

<Author's changes in manuscript>

l.1087-1100: "The second point is that the horizontal phase correlation has a large impact on nonharmonic internal-tide variance. As a simple example, consider the above two-source case but in the absence of horizontal correlation. Then, R^{1/2}=I and E(A'^2)=2|s_0|^2\varsigma_0^2 from Eq. (12), which is half of the above perfectly correlated cases. It is important to relate this to grid resolution in a numerical hydrodynamic model. If one source region is resolved by one grid point with s_phys=[2s_0] and \Sigma=\varsigma_0 in a low-resolution model and s_phys=[s_0 s_0]^T and \Sigma=\varsigma_0 I in the corresponding high-resolution model, the sum of s_phys (i.e., pre-modulation internal-tide amplitude) is the same (i.e., 2s_0). However, if we neglect the horizontal correlation of the sources, the variance is $E(A'^2)=4|s_0|^2\varsigma_0^2 in the low resolution case and 2|s_0|^2\varsigma_0^2 in the high resolution case. The perfect correlation considered in the last paragraph is required to make the variance the same at the two resolutions. This shows that the horizontal correlation has to be considered for gridded sources, otherwise the results would be highly dependent on grid resolution."

-------------------------

<Referee comments>

l253: Usually the "objective function" is a quadratic expression in data assimilation, so this is a little confusing. Why not just refer to J as an arbitrary linear function of x?

<Author's response>

Thank you for your comment. I followed your suggestion. I also deleted "objective function" as much as possible, except when they are related to data assimilation.

<Author's changes in manuscript>

-l.1120: "Using the model solution, we consider a linear function J=w^H x, tentatively defined at a particular time t_j."

-------------------------

<Referee comments>

Fig 2: Please put panel labels in the same relative location in each panel, i.e., at the top left.

<Author's response>

Thank you for your comment. The figures were modified following your suggestion.

<Author's changes in manuscript>

- Panel labels were moved to the top left in Figs. 3, 4, 7, 10 in the revised manuscript.

-------------------------

<Referee comments>

Fig 2d: Some of the structure depicted in this figure looks like it could be caused by spurious bottom topography data, e.g., such as the linear features around 119E, 19S and the apparent correlation of the forcing function with the 100m 200m and 500m isobaths. I would be curious to see a histogram of depths from this region to see if it exhibits peaks at 100m intervals.

<Author's response>

Thank you for your comment. I think it is worth pointing out that the data source for the Australian shelf in GEBCO 2019 bathymetry is the historical surveys compiled by Geoscience Australia (e.g., Whiteway 2009). From my experience in numerous commercial projects on the Australian North West Shelf, the Geoscience Australia bathymetry is a good choice as publicly available data, unless detailed surveys are available in the area of interest. This is one of the major reasons why GEBCO bathymetry is used in this study. (Updated products became available since I originally set-up the numerical model in 2019, but I believe the details are not important for this feasibility study.)

The reason why you see banded structure around 100m, 200m, and 500m isobaths is that the bottom slopes are steeper at these depth ranges in the region. Note that, from the definition in Eq. (A4), the forcing function is proportional to the vertically integrated barotropic transport (h_0 u_0, h_0 v_0) and the topographic interaction coefficients (L^x_{0n}, L^y_{0n}), which are proportional to bottom slope divided by water depth (see e.g., Eq. (B3) in Shimizu 2019). Also, the reason why you see apparent anti-correlation around (119E, 18.5S) is that the slope angle changes the sign because of a underwater ridge.

In the separately uploaded Fig. R1, I hope you see that the map of bottom slope divided by water depth has features similar to the forcing function (except the sign). Also, I hope you see that the bathymetry along the three transects show one "shelf break" less than 100m depth, and another shelf break or a ridge around 150m depth. They make the bottom slope steeper around 100m and 200m water depth. Also, the bottom slope tends to become steeper around 500 m, which increases bottom slope divided by water depth around 500m depth. (You probably notice fine orange lines in Fig. R1a. They are artefacts caused by imposing minimum bottom cell thickness for numerical modelling.)

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

l306: "same reasoning" refers simply to treating the finite-dimensional linear sum (eqn (12)) as an approximation of the integral? I am not convinced that the arguments about the matrix R translate directly to the function R, here.

<Author's response>

Thank you for your comment. You are right - I realized I forgot explaining one step. Because R is factorized as R^(1/2)R^(T/2), E(A'^2) can be written as the inner product of (R^(T/2)\Sigma s_phys) and its complex conjugate in the discretized form. So, the element-wise product can be used for mapping purposes. Therefore, the continuous version should be the integral of (integral (R^(1/2)\sigma s) times its complex conjugate).

I also forgot explaining that the nonharmonic variance source function is defined for the variance of time series for comparisons with observations and numerical modelling, instead of the variance of the envelope amplitude (the factor 1/2 comes from the variance of the "carrier" wave). The explanation was added in the revised manuscript.

This paragraph is moved to Section 3.1 (immediately after the discretized version, Eq. (12)), because the details of adjoint frequency response analysis was moved to new Appendix D.

<Author's changes in manuscript>

- The equation (Eq. (13) in the revised manuscript) was modified as explained above. The corresponding change was made to the discretized version, Eq. (12).

- The referred paragraph was moved to the end of Section 3.1, and the explanation of the equation was modified accordingly.

- l.231: "The factor 1/2 is multiplied in the above equation so that the integral of s_nh corresponds to the variance of nonharmonic internal-tide time series from observations or numerical modelling, rather than the variance of the envelope amplitude."

-------------------------

<Referee comments>

l310: Regarding the non-uniqueness: Aside from the explanation in 3.1, isn't it more fundamental that s_{nh} is not unique? s_{nh} is a function (it has inifintiely many degrees of freedom), while E(A'^2) is a scalar. Therefore it is not possible to uniquely determine s_{nh} from E(A'^2). This is distinct from the non-uniqueness of R^{1/2} discussed in 3.1.

<Author's response>

I am sorry that this was not clear. I meant to refer to only non-uniqueness resulting from the non-uniqueness of R^{1/2}. I do use tools developed for inverse modelling in this paper, but the model (with large degrees of freedom) was not constrained to the measurement (a scalar), as done in data assimilation.

I think that part of the problem was that I forgot explaining one step in the calculation of s_nh. With the revision in response to your last comment, I hope it is now clearer how the non-uniqueness of R^{1/2} affects s_nh.

<Author's changes in manuscript>

- Eq. (12) (discretized) and Eq. (13) in the revised manuscript (continuous) were modified as in my response to your last comment.

- l.236: "(However, note that s_nh is non-unique within the correlation length of phase modulation, because R^{T/2} in Eq. (12) or R^{1/2}(x',x) in Eq. (13) is non-unique, as explained in Appendix C.)"

-------------------------

<Referee comments>

Fig 3a: Does this figure contain a spurious pink line? There is a straight line at about 116.5E from 9S to 17S which looks out of place and does not appear to represent a ray path.

<Author's response>

I checked the results, and the lines are not spurious. (It happened that there were two paths that were very close. One of them was removed.) In the figure, the rays appear straight over Argo Abyssal Plain (~5700m water depth), because the abyssal plane is very flat and mean phase speed without background currents is used for ray tracing. Probably, you are more familiar with the results including phase-speed variation caused by mesoscale variability, such as Park and Watts (2006) and Rainville and Pinkel (2006). In response to your other comments, some justification for using ray tracing with mean stratification was added in Discussion.

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

Sect 3.3: While this section makes true statements, equation (26) is do general that I'm having trouble seeing how this will be used. Why not simply present (26) in the more specific context, where components of P, B, and Q are defined.

<Author's response>

Thank you for your comment and suggestion. After considering your and another referee's comments, I decided to move Section 3.3 in the original manuscript to new Appendix E, and the explicit forms of matrices A, B, and Q are shown for the two-path case immediately after the introduction of the covariance equations. Then, the equations used for the phase modelling were obtained by simplifying the covariance equations and the matrices.

<Author's changes in manuscript>

- Section 3.3 in the original manuscript was moved to new Appendix E.

- The explicit forms of matrices A, B, and Q are shown for the two-path case immediately after the introduction of the covariance equations (in Appendix E).

-------------------------

<Referee comments>

l350: Rather than refer to this as Lorentzian, which is usually applied to spectra of narrowband process with phase modulations, I wonder if it would be better to call it a first-order autoregressive process? I guess it is Lorentzian, but centered at the zero frequency.

<Author's response>

Thank you for your comment. Because the spectral form is not needed in this paper, I deleted the reference to Lorentzian spectrum.

<Author's changes in manuscript>

- The reference to Lorentzian spectrum was removed from l.1213 in the revised manuscript (now in Appendix E).

-------------------------

<Referee comments>

l354: P_{cc} is "stationary" -- do you mean P_{cc} constant in space and time?

<Author's response>

Yes, it needs to be stationary in space and time. The following is the equivalent sentence in the revised manuscript (in Appendix E).

<Author's changes in manuscript>

l.1205: "In this paper, we assume that the phase-speed variance P_{c_i c_j} is stationary in space and time as a first approximation (justified in Section 4.3)."

-------------------------

<Referee comments>

Eqns (31a-b): Are these derived from (27) and (28), which are the explicit form of equation (25)?

<Author's response>

I am sorry that this was not clear. Yes, (27) and (28) are put in the form of (25), and then the associated covariance equation (26) provides (31a-b). In response to your other comments, the explicit forms of matrices A, B, and Q are shown for the two-path case after the introduction of covariance equations (now in Appendix E). Then, the explicit form is simplified to derive Eq. (31a,b) (now Eq. (18)) as follows.

<Author's changes in manuscript>

l.1222: "Eq. (18) for the phase spread modelling is obtained from Eqs. (E2) and (E3) by neglecting the rows and columns corresponding to the ith path and the cross-path correlation (i.e., F = 0), and by writing \theta_j = \theta and c_j' = c'."

-------------------------

<Referee comments>

l363: The spatial variability of  \overline{c} and L_C is included in (31a-b), not in equation (32), right?

<Author's response>

Yes, that is correct. The sentence before Eq. (32) was revised as follows.

<Author's changes in manuscript>

l.344: "If \overline{c} and L_C remain constant, the solution under the initial condition P_{c\theta} = P_{\theta\theta} = 0 at t = 0 is ..."

-------------------------

<Referee comments>

l368: This is a little confusing. I think you are saying that for each source, j, there is an associated P_\theta\theta, the value of which is determined by the path between the source j and the observation location where P_\theta\theta (\sigma^j) is evaluated.

<Author's response>

I am sorry that the explanation was confusing. I think the problem was that the covariance equations provide the forward model, but my calculations were done backward in time, which was explained in Methods section. I moved sentences regarding the adjoint calculation of the covariance equations from Methods section to Theoretical background section. I also noted that P_{\theta\theta} grows with distance from the observation location.

<Author's changes in manuscript>

- l.352-364: "The straightforward approach for solving Eq. (18) is to integrate the equations from a source location to the observation location; however, this approach is computationally inefficient, because it needs separate (forward) integration from each source location along the same path. Alternatively, we can exploit the adjoint method described in Appendix D. The adjoint sensitivity of P_{\theta\theta} at the observation location to [P_{c\theta} P_{\theta\theta}]^T at other locations can be calculated by integrating the equations adjoint to Eq. (18) once, backwards in time from the observation location. Then, P_{\theta\theta} can be calculated as the convolution of the adjoint sensitivity and the forcing (i.e., \sigma_C^2 term in Eq. (18)) along the path. The resultant phase variance P_{\theta\theta}, which grows with distance from the observation location, is used as the phase variance \sigma_j^2 in the statistical model. Note that P_{\theta\theta} can grow without a limit, but this does not cause any problem because the wrapped normal distribution, Eq. (2), can be used with arbitrary large phase spread \sigma_j."

- l.406 "Similar to Eq. (18), Eq. (20) can be solved using the adjoint method explained in Appendix D, and the resultant variance P_{\Delta\theta\Delta\theta} corresponds to E(\Delta\Theta^2) in Eq. (7)."

- The following sentences were deleted.

-- l.367-370 in the original manuscript: "Note that Eq. (31) yields P_{\theta\theta} that increases from the source towards the observation location, but we associate the final P_{\theta\theta} with the source (the initial location of integration) in the statistical model. This is because we consider internal tides observed at a location in an "inverse" sense, and waves from remote sources are more random."

-- l.425-426 in the original manuscript: "(As in the case of P_{\theta\theta} with Eq. (31), the resultant P_{\Delta\theta\Delta\theta} is associated with the sources in the statistical model.)"

-------------------------

<Referee comments>

l377-l382: This is an interesting point. Presumably, though, there is an high-frequency cutoff. The whole ray-tracing and propagation paradigm only makes sense if c' is slowly varying compared to \omega. Right?

<Author's response>

Thank you for your comments. This is a detailed point but I think an important point to consider in the future. I do not have clear answers to the points you raise, because they appear open questions to me. However, some of my thoughts are provided below.

Regarding the high-frequency cutoff, I think there are several possibilities.

- The roll-off of power spectrum may be fast enough that high-frequency contributions may always be capped.

- In this study, the stochastic differential equation for c' (and the Lorentzian model in Part I) assumes only a single length scale for simplicity. However, in reality, there are processes with different length scales. Multiple scales may be accounted for by imposing high-frequency cut-off, or by using a model with multiple length scales. (Note that it is possible to introduce multiple scales in the stochastic equation for c'.)

Note that, in the analysis of the PIL200 observations (in Appendix A in the original Part I, which is now Appendix F in Part II), high-frequency cut-off at 1 hr period is effectively applied by low-pass filtering.

Regarding the slowly varying assumption, ray tracing formally requires the assumption. However, the answer to your question is not that simple, because what matters in the end is the phase statistics at the observation location, and their sensitivity to c' (ideally as a function of its length scale in a spectral sense). Also, I think there are many results we can learn from studies of "wave propagation in random media" in other fields of physics and engineering, such as acoustics (e.g., Colosi, 2016), regarding this point. For example,

- For a point source, the observed phase is insensitive to the small-scale variability of phase speed.

- Although c' can induce wide spread of internal-wave rays (Park and Watts, 2006; Rainville and Pinkel, 2006), the contributions to the observed statistics may come only from paths around the unperturbed propagation path (i.e., a Fresnel zone). This is because waves arriving through widely disturbed paths tend have different phases, and hence tend to average out through interference.

- Because of the above point, the unperturbed path may provide a convenient basis to calculate phase statistics.

Overall, I accept that there are many potential deficiencies with ray tracing, and I used ray tracing as a compromise. However, considering the above positive aspects, the approach used in this study (ray tracing using the mean stratification) does not appear unreasonable as a first approximation. I think the details need to be investigated in the future.

<Author's changes in manuscript>

Following paragraph was added in Discussion.

l.820-833: "Since the analysis in Appendix A suggests that nonlinear effects do not have the leading-order effects, the most important caveat of the proposed approach appears to be the use of ray tracing and the mean stratification to calculate wave propagation paths. The use of ray tracing may be questioned because, when phase-speed variability is included in ray tracing, the length scale of phase-speed variability can be comparable or shorter than the wavelength (invalidating the slowly varying assumption), and ray paths could vary widely (Park and Watts, 2006; Rainville and Pinkel, 2006). However, studies on wave propagation in random media in other fields, such as acoustics (e.g., Colosi, 2016), suggest that ray tracing may have wider applicability than it seems. For example, observed phase tends to be insensitive to small-scale phase-speed variability (consistent with Fig. 11a). Even when ray paths diverge widely, the contributions to the observed phase lag may come only from paths around the mean (unperturbed by phase-speed variability) propagation path, called a Fresnel zone. This is because waves arriving through widely perturbed paths tend to have different phases, and hence tend to average out through interference. They suggest that phase statistics has relatively weak dependence on the details of ray paths and small-scale phase-speed variability, which appears to be consistent with Buijsman et al. (2017). Ray tracing and the mean stratification are used in this study as a compromise among these factors and its simplicity. It would be worth investigating the impact of different methodologies for calculating wave propagation paths in the future."

-------------------------

<Referee comments>

Eqn (33): The last row of A, does it represent the evolution of \Delta\theta'', the phase evolution along the two paths? Is that why the terms are opposite signs, because \theta'' = \theta'(path 1) - \theta'(path 2)>

<Author's response>

Yes, it is. I am sorry that it was not clear. In response to your other comments, the explicit forms of matrices A, B, and Q are shown for the two-path case after the introduction of covariance equations (now in Appendix E). Then, the explicit form is simplified to derive Eq. (33) (now Eq. (20)) as follows.

<Author's changes in manuscript>

l:1224: "Eq. (20) for the cross-path phase difference modelling is obtained from Eqs. (E2) and (E3) by modifying the definition of x in Eq. (E1) as x = [c_i' c_j' theta_i-\theta_j]^T, and by subtracting the fourth row from the third row in A and B."

-------------------------

<Referee comments>

l419: The assumption of "time-independent" also mean "space-independent" in this context, becauase time is measured along ray paths. Is that correct?

<Author's response>

Yes, it is. Thank you for pointing it out. The following change was made in the revised manuscript.

<Author's changes in manuscript>

l.399: "However, if \overline{c}, L_C, and |\Delta\eta|/l remain constant, ..."

-------------------------

<Referee comments>

Eq (39): Please clarify: even though \overline{c}, L_C and \Delta\eta/l are assumed to be time-independent, P_\theta\theta is not time-independent and varies with t according to equation (32). But there is an oddity: the observation is at a signle point, and the path separation, \Delta\eta, presumably linearly decreases to 0 approaching this point. I wonder if there is a cancellation of the linear grown with time (equation (32)) and the linear decrease with time (\Delta\eta),.

Aha: l425-l440: It appears that you have already thought-through the consequences of my comment about Eq (39).

Overall comment on this section: I found the development a little hard to follow.

I wonder if it might be more direct to state the covariance evolution equations all at once, eqns (31)a-b and (38)a-c, and explain how these describe the along-path and across-path phase covariances (and c'). I'm not sure of the best approach. Perhaps it would be sufficient to add a short paragraph after eqn (26) with an overview of the approach to follow, so that the reader is prepared to accept the notation for along-path and acorss-path phase variations and their covariances. Maybe change the header of 3.4 to mention the "along-path phase difference", which would make it more parallel to the header of section 3.5.

<Author's response>

Thank you for your comments and suggestion. Considering your and another referee's comments, I have moved Section 3.3 in the original manuscript to new Appendix E, and the explicit forms of matrices A, B, and Q are shown for the two-path case (using x = [c_i' c_j' \theta_i' \theta_j']). Then, the covariance equations used in the stochastic phase modelling are introduced by simplifying this general case (by removing one of the paths for the phase spread modelling, and by changing x to x = [c_i' c_j' \theta_i'-\theta_j'] for the phase correlation modelling). The amount of text is about the same, but I hope this removed the need to consider the similar problems twice from the beginning.

<Author's changes in manuscript>

- Section 3.3 in the original manuscript was moved to new Appendix E.

- The explicit forms of matrices A, B, and Q are shown for the two-path case immediately after the introduction of the covariance equations.

- l.403: "(It may appear odd to assume constant |\Delta\eta|/l because |\Delta\eta| certainly vary; however, an empirical relationship is introduced later in Section 4.4 to account for the variation.)"

-------------------------

<Referee comments>

l449-l456: I'm sorry, but I don't understand what the phrase "were included as forcing of nonharmonic internal tides in the semidiurnal frequency band" means. Maybe you could start by clarifying what you mean by "in the modeling". Does this refer to ray tracing (modeling internal waves per se), or does it refer to some version of equation (26) (modeling phase/phase speed covariance), or some combination of these?

<Author's response>

I am sorry that this was not understandable. It is for the whole model suite. New Fig. 5 shows a flow chart for the application of the proposed model suite to multiple tidal constituents and vertical modes. Also, the paragraph was revised as follows. I hope it makes the paragraph understandable.

<Author's changes in manuscript>

-l.432-445: "In the model suite, we included the four major semidiurnal tidal constituents (M_2, S_2, K_2, and N_2) and four lowest baroclinic modes (VM1--VM4). Fig. 5 shows a flow chart for the application of the proposed model suite to multiple tidal constituents and vertical modes. Forcing from the major constituents were considered separately, assuming that the nonharmonic internal-tide variance (and the associated statistics) is calculated for a sufficiently long time series. Since it was impractical to separate nonharmonic internal tides into constituents in the PIL200 observations, the resultant variance, E(A'^2) in Eq. (13), and the nonharmonic variance source functions from individual constituents were summed to obtain the total for semidiurnal internal tides. It may sound confusing to include multiple baroclinic modes to model VM1 internal tides at the PIL200 location. This is required because barotropic forcing excites not only VM1 but also higher modes, which can be converted to VM1 by topographic interaction before arriving at the PIL200 location (see Fig. 5). To distinguish overall barotropic forcing to VM1 internal tides at the PIL200 location from barotropic forcing to individual baroclinic modes in the intermediate process, the latter is hereafter referred to as, for example, "barotropic-to-VM2" or "VM0-to-VM2" forcing."

-------------------------

<Referee comments>

l460: omit "in this feasibility study" --- We are 18 pages into this manuscript, and a number of new concepts have been introduced. You need to simplify the language as much as possible to make this story comprehensible.

<Author's response>

Thank you for the suggestion. I deleted "in this feasibility study" from Section 2.

<Author's changes in manuscript>

- l.460 in the original manuscript: "in this feasibility study" was deleted.

- l.498 in the original manuscript: "In this feasibility study" was deleted.

-------------------------

<Referee comments>

l464: "amplitude" --- Need to be careful here. The amplitude, sqrt(A^2 + B^2), is not a linear functional of the time series A cos(omega t) + B sin(omega t). I'm not sure if the model you are describing is for the waves, or for their statistics. If it is the latter, then the amplitude could very well be a linear functional of the state.

<Author's response>

Thank you for your comment. I meant vertical-mode amplitude, or the complex-valued amplitude of the time series, a e^(omega t-\phi). I added the word "complex-valued" and its mathematical expression "a e^{-i\varphi}" in the revised text. Also, it appears that this comment is related to another comment regarding the use of the term "objective function". I removed the phrase "objective function J" from the sentence.

<Author's changes in manuscript>

l.455: "It was calculated for complex-valued VM1 isopycnal-displacement amplitude at the PIL200 location (i.e., a e^{-i\varphi} in Eq. (14)), whose magnitude was scaled to have the value of extreme (maximum or minimum) displacement within the water column."

-------------------------

<Referee comments>

l470: The "model" here clearly describes the model for the wave propagation.

<Author's response>

Thank you for your comment. I added "hydrodynamic" in front of "model" in two sentences in the subsection.

<Author's changes in manuscript>

- l.453: "A sinusoidal periodic motion was assumed (as in Eq. (D7)) in Appendix D) in the governing equations (Eq. (B3) in Appendix B without the nonlinear terms), so that the hydrodynamic model directly calculates the adjoint frequency response function (\lambda in Eq. (14))."

- l.462: "Details of the hydrodynamic model set-up were as follows."

-------------------------

<Referee comments>

l479-l480: I don't understand what this means. How does the assumption that the barotropic and VM1 kinetic energies are the same relate to bottom friction? Or, is all of l479-l484 to be interpreted together (but this seems to contradict the simpler explanation in l478)?

<Author's response>

I am sorry that this was not clear. Because this is a detail of the model set-up and bottom friction has substantial effects only on relatively shallow shelves, I retained only the simpler explanation in the revised manuscript, to make the manuscript shorter.

<Author's changes in manuscript>

l.477-482 in the original manuscript were deleted.

-------------------------

<Referee comments>

l495: Using the language of Bennett, it sounds like you computed representer functions and their adjoints (which you referred to as the sensitivities) for VM1-4 and the 4 tidal frequencies, assuming a measurement at PIL200.

<Author's response>

With reference to Chapter 1.3.3 in Bennett (2002) and Ebgert and Evofeeva (2002), I believe "source functions" in the manuscript are not representer functions, because I did not use the adjoint solutions to force the forward model. The adjoint solutions may be referred to as the adjoint of representer functions in Bennett's terminology, but it appears more standard to call them the Green's functions (Chapter 1.1.4), because I used harmonically analyzed point measurement (i.e., "initial condition" of the adjoint model is delta function). The source functions are the products of the respective Green's functions and internal-tide generation forces (i.e., forcing functions).

<Author's changes in manuscript>

No change was made to the manuscript based on this comment.

-------------------------

<Referee comments>

l502-l508: I would suggest moving this to the discussion.

<Author's response>

Thank you for your suggestion. Considering your other comments, I added a paragraph in Discussion regarding potential caveats of ray tracing. Since the new paragraph is not short, I decided to delete the material in l.502-508 in the original manuscript.

<Author's changes in manuscript>

The following sentences were deleted.

l.502-508 in the original manuscript: "This ray tracing method had potential deficiencies, such as the neglect of interaction of vertical modes, wavelengths that are not much shorter than the continental slopes, the existence of multiple paths to some regions, and the difficulty in calculating paths passing through straits. However, these potential deficiencies were considered to be relatively minor to the overall results of this study. The reason is that travel time and the phase variance from Eq. (31) have relatively weak dependence on the details of ray paths, because they are integrated quantities of spatially variable time-mean variables, such as phase speed and its correlation length."

-------------------------

<Referee comments>

l509-l516: How were \sigma_C^2 and L_c were chosen to integrate (31)? Aha -- I see later.

<Author's response>

Thank you for your comment. I think the problem was that there was another paragraph before the explanation of \sigma_C^2 and L_c. The referred paragraph was deleted to make the manuscript shorter, so I hope you do not have the same problem.

<Author's changes in manuscript>

The referred paragraph (l.509-l.516 in the original manuscript) was deleted.

-------------------------

<Referee comments>

l522: Are you saying that the non-M2 tides will use the same along-path (P_\theta\theta, P_cc, and P_c\theta) statistics, as determined for M2? You haven't mentioned yet the across-path statistics, but I assume those a included also?

<Author's response>

Yes. I have tried calculating phase statistics using the S2, K2, and N2 frequencies, but the results were essentially the same as those with the M2 frequency. It was the same for cross-path statistics. The sentence was modified as follows.

<Author's changes in manuscript>

l.515: "Since the results were insensitive to small frequency differences among the major semidiurnal constituents, the M_2 frequency was used in the modelling."

-------------------------

<Referee comments>

l551: "length scale larger than eddies" --- Indeed. You might want to cite Buijsman's study of the near-equator variability in HYCOM.

<Author's response>

Thank you for your suggestion. I added reference to Buijsman et al. (2017).

<Author's changes in manuscript>

l.545: "However, phase-speed correlation could be affected by processes that have length scale larger than eddies (e.g., Buijsman et al.,2017)."

-------------------------

<Referee comments>

l588-l589: I don't understand the parenthetical comment about "one dimensional".

<Author's response>

Thank you for your comment. Since it is unnecessary, the parenthetical comment was deleted.

<Author's changes in manuscript>

l.586: "For example, the correlation function is highly anisotropic for \sigma_\xi = 9\sigma_\eta, but it yields \alpha_r = 3."

-------------------------

<Referee comments>

l592-594: Why is this other estimate "more accurate"? Simplify if possible. And don't consider all the caveats (or move them to the Discussion).

<Author's response>

Thank you for your comment. Since it is not important, the sentences were deleted.

<Author's changes in manuscript>

The following sentences were deleted.

l.592-594 in the original manuscript: "A more accurate way to evaluate $\alpha_r$ is to consider the path average of the Gaussian correlation function calculated tentatively with \alpha_r=1, and then determine the equivalent \alpha_r for the same \Delta r. This shows that \alpha_r=2 is a reasonable choice."

-------------------------

<Referee comments>

l597-l608: I don't understand the significance of all these details. Perhaps this could be expanded somewhat to explain, or maybe all this could be pushed into the discussion or an appendix.

<Author's response>

I am sorry that this section was not understandable. I modified the paragraph substantially, and added steps taken to calculate the correlation length \sigma_r. I hope the details are understandable now.

<Author's changes in manuscript>

The paragraph was revised as follows.

l.595-611: "Based on the above consideration, the equivalent isotropic correlation length of phase modulation \sigma_r was calculated from Eqs. (7), (20), and (26) as follows. Considering both anisotropy of the phase correlation and the along-path variation of cross-path distance, \alpha_r between 1 and 5 appeared to be reasonable. We arbitrarily chose $\alpha_r=3$ as a reference value. Using Eq. (26) with a chosen \alpha_r, P_{\Delta\theta\Delta\theta} from Eq. (20) was substituted into Eq. (7) to calculate the correlation coefficient for different source distance \Delta r. This yielded isotropic correlation function R(\Delta r). Since \sigma_r was required for the diffusion operator method, the Gaussian shape was fit to the first peak of the correlation function where R(\Delta r) > 0.5 by the least-squares method, and the resultant standard deviation was used as \sigma_r in the diffusion operator method."

-------------------------

<Referee comments>

l616-l625: Hmmm. It seems like you should have sorted out the significance of these processes with some wave model runs using variable phase speed modulations. It sounds like you have made some arbitrary assumptions about the nature of the phase variability.

<Author's response>

I am sorry that this section was not understandable. In response to your another comment, new Fig. 5 was added to show a flow chart for the application of the proposed model suite to multiple tidal constituents and vertical modes. I hope this section is understandable with the figure.

I did make an assumption that higher modes are converted to VM1 near the PIL200 location, but in this particular example application, this is much more certain than the three model parameters. The reason is that the adjoint frequency functions for higher modes (e.g., Fig. 5b in the original manuscript) show dominant topographic conversion from VM1 to higher modes near the PIL200 location (note that the process is in the opposite direction because the adjoint model runs backwards in time). I think it is important to mention this assumption because it may not hold in other applications.

<Author's changes in manuscript>

- Reference to new Fig. 5 was added in the paragraph.

- l.630: "The latter scenario was assumed in this study, because the continental slope near the PIL200 location induced strong topographic interaction between VM1 and higher modes, as shown later."

-------------------------

<Referee comments>

l660: Remind us what is the "source function". Which term in which equation is it?

<Author's response>

Thank you for your comment. I added reference to the definition when adjoint frequency response function and source function are referred to in the Methods and Results sections.

<Author's changes in manuscript>

- l.453: "A sinusoidal periodic motion was assumed in the governing equations (as in Eq. (D7) in Appendix D), so that the hydrodynamic model directly calculates the adjoint frequency response function (\lambda in Eq. (14))."

- l.486: "The source function (s(x) in Eq. (14)) was calculated from ..."

- l.648: "The adjoint frequency response function (\lambda in Eq. (14)) of VM1-induced isopycnal displacement at the PIL200 location to the barotropic (VM0)-to-VM1 forcing qualitatively shows ..."

- l.667: "The source function (s(x) in Eq.14)) was calculated simply by ..."

-------------------------

<Referee comments>

Fig 8: Sorry if I have fogotten: why is it necessary to make a Gaussian fit to the covariance functions? Is it because the Weaver diffusion model is used with the wave evolution equation model? Or, do you just need to estimate the correlation length?

<Author's response>

Because the correlation length (standard deviation of Gaussian function) is required in Weaver and Courtier's diffusion operator method. I am sorry this was unclear. I think this comment is related to your earlier comment that the details regarding the calculation of the correlation length \sigma_r were not understandable. I added a reminder at the beginning of Section 5.4, and also in the caption.

<Author's changes in manuscript>

- l.688: "Since the diffusion operator method by Weaver and Courtier (2001) was used to represent the horizontal correlation of phase modulation, the equivalent isotropic correlation length of phase modulation \sigma_r characterises the horizontal correlation."

- Caption of Fig. 9: "Dotted vertical lines indicate standard deviations determined by least-squares fit of Gaussian function, which is used as correlation length \sigma_r for diffusion operator method by Weaver and Courtier (2001)."

-------------------------

<Referee comments>

l783-l785: I don't understand the distinction between the "total modelled variance" and the "VM1 M2 component". Does the total mean that all VM1-4 and the 5 harmonic constants are included? The amplitudes of the VM2-4 and non-M2 are so much smaller than VM1 M2, though, I'm not sure of the usefulness of the degrees of freedom concept for these smaller components.

<Author's response>

Thank you for your comment. Considering overall comments from you and another referee, I decided to remove all materials related to the PDF comparison and DoF to make the manuscript shorter. So, these sentences were deleted.

<Author's changes in manuscript>

The following material related to the PDF comparison and DoF were deleted.

- l.128 in the original manuscript: "The method for calculating probability density function (PDF), which is used only briefly near the end of this paper, is described in Part I."

- Section 4.6.

- The last paragraph in Section 5.6.

- The DoF columns in Table 1.

- Panel (c) and (d) in Fig. 10 in the original manuscript (now Fig. 11).

-------------------------

<Referee comments>

l821: Didn't Rainville and Pinkel's study find that the variability of propagation paths is "large", in some sense. At that time, there was some discussion of how this invalidated some aspects of the ray-tracing approach, but I cannot reconstruct the arguments. In any case, you proposed model (without the path variability) yields plausible results. But maybe this is just another aspect of the non-observabilkity of the detailed mechanisms in this problem.

<Author's response>

[The following response is mostly duplicate of my response to your earlier comment. But because there are many comments, I am repeating my response.]

I am aware that the variability of ray paths is large when phase-speed variability is included in ray tracing. However, the answer to your question is not that simple, because what matters in the end is the phase statistics at the observation location, and their sensitivity to c' (ideally as a function of its length scale in a spectral sense). Also, I think there are many results we can learn from studies of "wave propagation in random media" in other fields of physics and engineering, such as acoustics (e.g., Colosi, 2016), regarding this point. For example,

- For a point source, the observed phase is insensitive to the small-scale variability of phase speed.

- Although c' can induce wide spread of internal-wave rays (Park and Watts, 2006; Rainville and Pinkel, 2006), the contributions to the observed statistics may come only from paths around the unperturbed propagation path (i.e., a Fresnel zone). This is because waves arriving through widely disturbed paths tend have different phases, and hence tend to average out through interference.

- Because of the above points, the unperturbed path can provide a convenient basis to calculate travel time or phase statistics.

Overall, I accept that there are many potential deficiencies with ray tracing, and I used ray tracing as a compromise. However, considering the above positive aspects, the approach used in this study (ray tracing using the mean stratification) does not appear unreasonable as a first approximation. I think the details need to be investigated in the future.

<Author's changes in manuscript>

Following paragraph was added in discussion.

l.820-833: "Since the analysis in Appendix A suggests that nonlinear effects do not have the leading-order effects, the most important caveat of the proposed approach appears to be the use of ray tracing and the mean stratification to calculate wave propagation paths. The use of ray tracing may be questioned because, when phase-speed variability is included in ray tracing, the length scale of phase-speed variability can be comparable or shorter than the wavelength (invalidating the slowly varying assumption), and ray paths could vary widely (Park and Watts, 2006; Rainville and Pinkel, 2006). However, studies on wave propagation in random media in other fields, such as acoustics (e.g., Colosi, 2016), suggest that ray tracing may have wider applicability than it seems. For example, observed phase tends to be insensitive to small-scale phase-speed variability (consistent with Fig. 11a). Even when ray paths diverge widely, the contributions to the observed phase lag may come only from paths around the mean (unperturbed by phase-speed variability) propagation path, called a Fresnel zone. This is because waves arriving through widely perturbed paths tend to have different phases, and hence tend to average out through interference. They suggest that phase statistics has relatively weak dependence on the details of ray paths and small-scale phase-speed variability, which appears to be consistent with} Buijsman et al. (2017). Ray tracing and the mean stratification are used in this study as a compromise among these factors and its simplicity. It would be worth investigating the impact of different methodologies for calculating wave propagation paths in the future."

-------------------------

<Referee comments>

l922: Do you know if this system is equivalent to Sam Kelly's Coupled Shallow Water (CSW) representation? It looks equivalent except he used z (rather than sigma) as the vertical coordinate. Anyway, it might be worth mentioning the equivalence. Also, his formulation admits somewhat simpler expressions for the coupling coefficients (L_{nm}) than appears to be the case here (documented in an appendix of Zaron, Musgrave, and Egbert). Some earlier work by Lahaye has similar expressions for the mode-coupling terms. 

<Author's response>

For linear hydrostatic problems (including this study), the system is equivalent to Kelly's Coupled Shallow Water (CSW) model (at least in its original version in 2016), because the formulation of CSW and that in the appendix are based on Shimizu (2011), although Kelly and coworkers stopped citing the paper after Kelly et al. (2011). Zaron and Ebgert (2014) also used equivalent formulation. The simpler form of topographic interaction coefficients (L_{nm} or T_{nm}) in Zaron, Musgrave, and Egbert (2022) is equivalent to Eq. (B5) in Shimizu (2011), or its continuous version, Eq. (B3) in Shimizu (2019). The major difference between different coordinate systems appears in nonlinear problems. The isopycnal coordinate (Shimizu 2017, 2019) allows the extension of the formulation to full nonlinearity, whereas it is difficult to avoid weak nonlinearity and/or slowly varying assumption in the height (z) coordinate (e.g., Kelly et al. 2016).

Note that, although Shimizu (2011) uses a multi-layer model instead of continuous stratification, the difference is not fundamental because the resultant modal evolutionary equations are the same, as pointed out in the paper. Although I now see some shortcomings of the paper due to the lack of experience back then, the paper is under-rated because the ideas proposed in the paper (the evolutionary equations of modal amplitudes in a form analogous to the shallow water equations, the associated modal energetics in a conservative form, the explicit expression of topographic interaction coefficients, etc.) have been used in many studies after 2011, although most of these studies did not cite my paper.

If you compare the formulation in Shimizu (2011) to the earlier formulation by Griffith and Grimshaw (2007), the evolutionary equations are clearly different, although the two versions should be equivalent. For example, Eq. (11) in Griffith and Grimshaw (2007) is integro-differential equations including second-order derivatives of the modal structure function in the interaction coefficients, whereas the evolutionary equations in Shimizu (2011) are differential equations analogous to the shallow water equations, including only first-order derivatives in the interaction coefficients. The simpler formulation appears to have contributed to wider use of the method, because studies after 2011 adopted my formulation.

<Author's changes in manuscript>

l.1025: "The linear formulation by Shimizu (2011) was adopted by, for example, Zaron and Egbert (2014) and Kelly et al. (2016)."

-------------------------

<Referee comments>

l37: "tidal currents" --> "barotropic tidal currents"

l47: " nonstationary" --> "nonstationary"

l655: "sound nothing" --> "sound like nothing"; but you should probably rephrase this in the positive: "This observation is significant because ..."

<Author's response>

Thank you for finding errors and your suggestion.

<Author's changes in manuscript>

Errors were fixed, and the expression was modified following the referee's suggestion. 

-------------------------

<Referee comments>

l705-l728: This is nice qualitative discussion. Also, the following discussion of phase correlations is good, too.

<Author's response>

Thank you.

-------------------------

References

Kelly, S. M., Lermusiaux, P. F. J., Duda, T. F., and Haley, P. J.: A Coupled-Mode Shallow-Water Model for Tidal Analysis: Internal Tide Reflection and Refraction by the Gulf Stream, J. Phys. Oceanogr., 46, 3661–3679, doi:10.1175/JPO-D-16-0018.1, 2016. 

Kelly, S. M., Nash, J. D., Martini, K. I., Alford, M. H., and Kunze, E.: The cascade of tidal energy from low to high modes on a continental slope, J. Phys. Oceanogr., 42, 1217–1232, 2012.

Whiteway, T.: Australian Bathymetry and Topography Grid, June 2009. Scale 1:5000000. Geoscience Australia, Canberra. http://dx.doi.org/10.4225/25/53D99B6581B9A, 2009.

Zaron, E. D., Musgrave, R. C., and Egbert, G. D.: Baroclinic tidal energetics inferred from satellite altimetry, J. Phys. Oceanogr., 52, 1015-1032, doi: 10.1175/JPO-D-21-0096.1, 2022.