Wave propagation with irregular dissipation and applications to acoustic problems and shallow waters

In this paper we consider an acoustic problem of wave propagation through a discontinuous medium. The problem is reduced to the dissipative wave equation with distributional dissipation. We show that this problem has a so-called very weak solution, we analyse its properties and illustrate the theoretical results through some numerical simulations by approximating the solutions to the full dissipative model for a particular synthetic piecewise continuous medium. In particular, we discover numerically a very interesting phenomenon of the appearance of a new wave at the singular point. For the acoustic problem this can be interpreted as an echo effect at the discontinuity interface of the medium.


Introduction
This work is devoted to the investigation of the 1D wave propagation through a medium with positive piecewise regular density and wave speed functions. For these non-smooth data we show that the problem has a so-called 'very weak' solution. This notion has been introduced in [GR15] in the analysis of second order hyperbolic equations, and in [RT16b] it was applied to show the well-posedness of the wave equations for the Landau Hamiltonian with irregular electro-magnetic fields. In this paper we use it to prove the well-posedness of the acoustic problem. Moreover, it allows us to derive the decay properties in time also in the situation when the medium has discontinuities. Incidentally, the same model equation (see (1.4)) appears also in the shallow water equations (see (1.6)) as a special case of the linear Boussinesq system, so the obtained results apply in that situation as well.
We start with a description of the physical problem. Now, let t denote the time and let z be the Cartesian coordinate in the direction of the wave propagation. Let ρ denote the density and c the wave speed of the medium. Following the derivation in [BB95], we obtain the first order hyperbolic system with the constitutive equation and the momentum equation: (1.1) p t + ρc 2 ω z = 0, ρω t + p z = 0, where p = p(z, t) is the z-component of traction across surfaces z = C (C -positive in compression), and ω = ω(z, t) is the z-component of particle velocity at the point (z, t). The first equation in (1.1) is Hooke's law differentiated in t. By denoting we rewrite (1.1) with respect to the impedance ζ(x) = ρ(z)c(z) (see [BB83] for the detailed argument) and obtain: (1.2) p t + ζω x = 0, ζω t + p x = 0.
Putting together equations (1.2), we get the equation Putting the initial conditions at x = 0, the problem (1.2) was analysed in [BB95] giving a rigorous estimate of the error occuring in making the O'Doherty-Anstey approximation, originally derived in 1971 in the context of acoustic wave propagation in the earth's crust [ODA71]. Furthermore, in [BB95] it was shown that the downgoing wave component D = ζ 1/2 ω + ζ −1/2 p decays as x → ∞ (i.e. large propagation distance) for smooth positive functions ρ and c, so that ζ is also smooth and positive. An extensive physical discussion of this kind of model equations was done in [BB95], [BB83], [CF94], [LB96], [PS00], [AKPP91]. Extension of the O'Doherty-Anstey approximation for weakly-dispersive, weakly nonlinear water waves propagating on the surface of a shallow channel with a random depth was developed in [MN04], [GMN07], and [MN05].
In this paper we are interested in the problem of existence of solutions of the model equation (1.3) in the situation when the density ρ and the wave speed c of the medium are irregular. For example, we want to allow them to be discontinuous or, in general, have even less regularity. At the same time, it is natural, from the physical meaning of these functions, to continue assuming that they are positive: ρ > 0, c > 0 so that also ζ = ρc > 0.
The first idea in our analysis is to observe that if we change the roles of t and x, the equation (1.3) takes the form of the dissipative wave equation. Indeed, swapping variables x and t, and denoting u := ω and b := ζ, the boundary value problem (1.3) is reduced to the Cauchy problem for the one dimensional acoustic equation has to be positive if we want to have the decay of solutions for large t. Indeed, this corresponds to the system losing energy (rather than gaining it from outside) and is a natural physical assumption. Since b > 0, this means that we should also have b ′ > 0. From the theory of distributions we know that these conditions would imply that b and b ′ are positive Radon measures, so that it is natural to assume that b is a piecewise continuous (and hence also increasing) function. Therefore, these will be the assumptions for our analysis, namely, we assume that the product b = ζ = ρc of the density and the wave speed of the medium is an increasing piecewise continuous function. If they are smooth (or at least C 1 ), this is a natural physical assumption (see [BB95]). Thus, the main novelty of this paper is that we relax the regularity assumption (1.5) requiring only that b is piecewise continuous, positive and increasing.
While physically this is a very reasonable setting, mathematically we face several problems: is not well-defined as a distribution; • moreover, if the data u 0 and u 1 are irregular, the dissipation term b ′ (t) b(t) ∂ t u(t, x) also does not make sense as a distribution in view of the celebrated Schwartz' impossibility result in [Sch54] on multiplication of distributions.
Nevertheless, in this paper we analyse the equation (1.4) or, more generally (2.1), under the assumption (1.5). In particular, we show that it is well-posed in the sense of very weak solutions introduced in [GR15], and then also used in [RT16b] in another context. In particular, we show that • if the Cauchy data (u 0 , u 1 ) is in the Sobolev spaces H s+1 × H s , s ≥ 0, then the Cauchy problem (1.4) has a very weak solution of order s; the very weak solution is unique in an appropriate sense; • if b ∈ C 1 we know that the Cauchy problem (1.4) also has a classical solution in C([0, ∞), H s+1 )∩C 1 ([0, ∞), H s ); in this case the very weak solution recaptures the classical solution; • we also give the above results for the Cauchy data (u 0 , u 1 ) in the Sobolev spaces H s+1 × H s for any s ∈ R; • under assumption (1.5), the very weak solution is uniformly bounded in L 2 and may decay in t depending on further properties of b.
The last property of the decay depending on further properties of b is natural even for smooth functions b, see [Wir04,Wir06,Wir07]. We note that if b is regular, we could write b ′ (t) b(t) = (log b(t)) ′ , which could be interpreted as a measure since log b(t) is a function of bounded variation. However, such a relation appears to be only formal since the quotient b ′ (t) b(t) is not well-defined even as a distribution. This also gives a novelty compared to the setting of [GR15] since the coefficients of the wave equations there were assumed to be distributions. We also note that for functions b more regular than C 1 , there are many results available, for the well-posedness in Gevrey spaces and in spaces of ultradistributions, see e.g. [CDGS79,DS98,KS06] and references therein, to mention only very few.
Let us mention that the equation (1.4) also appears in the modelling of the nondispersive water wave propagation in shallow water. Namely, the linear case of the Boussinesq systems takes the form Indeed, this is the non-dispersive case of the system analysed in [MN04, (2.5)] (the case of β = 0) to which we refer to further physical details. The function M is determined by the channel depth and local wave speed, and is discontinuous in channels with sudden changes in the depth.
Again, with the change of variable where C 0 (x) = 1/M(x), we arrive at the equation which is the same model (1.3), where the coefficient C 0 (the local wave speed) has the role of the impedance.
Numerics. In Section 6 we make a numerical modelling of the problem with b(t) having a jump discontinuity. We see that the approximating technique of very weak solutions does allow us to recover some physically expected behaviour such as the decay of solutions for large times, making such considerations mathematically rigorous. Moreover, we observe a very interesting phenomenon of the appearance of a new wave after the singular time travelling in the direction opposite to the main one. For the original acoustic problem this can be interpreted as an echo effect at the discontinuities of the medium. Such phenomenon is known in the presence of multiple characteristics in hyperbolic equations and is related to conical refraction, see e.g. [MU79] or [Lie93] for different descriptions. In that case the newly appearing wave is weaker, see e.g. [KR07], which is consistent with the observed pictures. The difference, however, is that in our case such behaviour is not related to multiple characteristics but to the singularity in the dissipation coefficient. Indeed, we analyse it further in Figures 5-7 by modelling a delta distribution as the Cauchy data. We see that the second wave is not more regular than the main one: it appears to be the second delta function (thus, singularity of comparable strength to the first one) but smaller in amplitude.
In Section 2 we discuss the known results in the case when b is regular, and then formulate our results concerning very weak solutions for irregular b. We also consider the problem in higher dimensions, which is of independent mathematical interest. In Section 4 we prove the results. In Section 5 we briefly discuss the case of distributional Cauchy data. Finally, in Section 6 we include some numerical experiments to illustrate the theoretical results derived in this paper for a particular synthetic piecewise continuous impedance function ζ.

Main Results
In this work we can deal with the Cauchy problem (1.4) in R n since our approach works in any dimension n ≥ 1. More precisely, we can consider where ∆ is the Laplace operator in R n . From the point of view of the acoustic problem the case n = 1 is relevant, but mathematically the equation (2.1) with irregular b as in (1.5) is of independent interest. We start by briefly recalling the notion of very weak solutions. For this, we do not need the assumption on the continuity of b since it is available in a much more general context.
As already mentioned, we will be using the notion of very weak solutions which was formulated for wave equations with spacially constant coefficients in [GR15], and applied to the wave equation for the Landau Hamiltonian in [RT16b]. We start by regularising the distributional coefficient b with a suitable mollifier ψ generating families of smooth functions (b ε ) ε , namely and ω(ε) is a positive function converging to 0 as ε → 0 to be chosen later (sometimes we need a particular behaviour in ε, see [GR15] and [RT16b] but the situation for (2.1) is simpler). Here ψ is a Friedrichs-mollifier, i.e. ψ ∈ C ∞ 0 (R), ψ ≥ 0 and ψ = 1. It follows that the net (b ε ) ε is C ∞ -moderate, in the sense that its C ∞ -seminorms can be estimated by a negative power of ε. More precisely, following [GR15], we will use the notions of moderateness: Here and in the sequel, the notation K ⋐ R means that K is a compact set in R. We note that the conditions of moderateness are natural in the sense that regularisations of distributions are moderate, namely we can regard by the structure theorems for distributions. Following and adapting [GR15], we now define the notion of a 'very weak solution' for the Cauchy problem (2.1): In the sequel, the proofs of the main statements (and especially the decay of the very weak solutions) will depend on behaviour In order to put the problem in perspective and for our subsequent use, we consider several cases. First, when for a fixed ε > 0, we have t(b ′ ε (t)/b ε (t)) → ∞ as t → ∞ then we will use the following result following from [Wir07]: ) → ∞ as t → ∞, then the solution to (2.4) and its derivatives satisfy the Matsumura-type estimates and as t → ∞ for any s > 0 and for all α ∈ N n 0 . Theorem 2.3 follows from [Wir07, Result 1] by setting p = q = 2. Integrals in (2.5) and (2.6) are well-defined since b ′ ε is also a positive function which will be discussed below.
In the cases when t as t → ∞ we will use the following results from [Wir06]: Then the solution to (2.4) satisfies the estimates and as t → ∞ for arbitrary s > 0.
We will also need the following extension of Theorem 2.4, keeping its assumptions: Corollary 2.5. The solution to (2.4) satisfies the estimates as t → ∞, for arbitrary α ∈ N n 0 for all s > 0 and for l = 0, 1.
Then for these functions we consider the Cauchy problem (2.10) Therefore, by Theorem 2.4 the solution v ε of the equation (2.10) satisfies and and Proposition 2.6. Let b be a positive distribution and let ε > 0. Consider the following two cases: Then, respectively, for all s > 0 we have (a) the solution of the Cauchy problem (2.4) satisfies the estimate for some constant C which is not depending on b, ε and initial data; (b) the solution of the Cauchy problem (2.4) has the decay as t → ∞, for some constant C which is not depending on b, ε and initial data. If b is a positive, piecewise continuous and increasing function, then the solution of the Cauchy problem (2.4) is uniformly bounded in ε ∈ (0, 1], i.e. for the solution u ε the following estimate is true where the constant C b,u 0 ,u 1 depends only on b and initial data u 0 , u 1 .
Proof. Lemma 2.7. Let us consider the nonhomogeneous equation Then, for any t ∈ [0, T ] and s ∈ R we have Proof. We begin by applying the Fourier transform with respect to the spatial variables. This yields the Cauchy problem for the second ordinary differential equation (2.14) Now, instead of (2.14) consider the following system with the Cauchy data The matrix (2.17) is symmetric. This allows us to use statements of e.g. Taylor's book [Tay81, Chapter IV] (also see [GR15b, Case 1] and [GR12]). Hence, we obtain the statement of the lemma. Now let us formulate the main results of this paper. As described in the introduction it will be natural to make the assumptions (1.5) for our analysis, ensuring that b ′ /b is positive. Theorem 2.9 (Consistency). Assume that b ∈ C 1 ([0, ∞)) is an increasing function such that b ≥ b 0 > 0, and that b ′ is a positive function such that b ′ ≥ b ′ 0 for some constant b ′ 0 > 0. Let s > 0, and consider the Cauchy problem with (u 0 , u 1 ) ∈ H s × H s−1 . Let u be a very weak solution of (2.18). Then for any regularising family b ε in Definition 2. We note that the convergence in Theorem 2.9 can be realised also on the interval [0, ∞) depending on further properties of b(t) allowing a global application of Lemma 2.7, see Remark 4.3.

Proof of Theorem 2.8
We regularise b by the convolution with a mollifier in C ∞ 0 (R) and get nets of smooth functions as coefficients. More precisely, let ψ ∈ C ∞ 0 (R), ψ ≥ 0 with ψ = 1. Define By the structure theorem for compactly supported distributions, we have that there exist L ∈ N 0 and c > 0 such that for all k ∈ N 0 and t ∈ [0, T ]. We note that the numbers L may be related to the distributional orders of b.
Now, let us define L more precisely under the assumption (1.5). Indeed, when b is a piecewise continuous, positive and increasing function, we obtain for all t ∈ [0, T ], where c, c 1 depend only on T . Thus, for b, which is a piecewise continuous, positive and increasing function, L = 0, and for the distributional function (δ-like) b ′ , we have L = 1.
Now the proof is depending on the behavior of b ′ ε (t)/b ε (t) at t → ∞. Let us consider several cases: when tb ′ ε (t)/b ε (t) → ∞ as t → ∞, and when tb ′ Anyway applying Theorem 2.3, or Theorem 2.4, or Corollary 2.5 to the equation (3.2), using the inequality (3.1) and that 1 |b we get the estimate The fact that there exist N ∈ N 0 , c > 0 and, for all k ∈ N 0 there exist c k > 0 such that ∂ k t u ε (t, ·) H s−k ≤ c k ε −N −k , for all t ∈ [0, ∞), and ε ∈ (0, 1], follows from the estimate which holds for all t ∈ [0, ∞), k ∈ N 0 , and ε ∈ (0, 1], from Theorems 2.3 and 2.4, Corollary 2.5, and acting by the iterations of ∂ t and by ∆ x on the equality and taking it in L 2 -norms. It means that u ε is from the space This shows that the Cauchy problem (2.1) has a very weak solution.

Proof of Theorem 2.9
Here we prove the consistency of the very weak solution with the classical one when the coefficients are regular enough. But first, we show that the very weak solution is unique in an appropriate sense. To present uniqueness we will use the notions of Colombeau algebras adapted to the properties of classical solutions.
Definition 4.1. We say that (u ε ) ε is C ∞ -negligible if for all K ⋐ R, for all α ∈ N and for all ℓ ∈ N there exists a constant c > 0 such that for all ε ∈ (0, 1]. Also, we call (f ε ) ε is C ∞ (R + ; H s )-negligible if for any K ⋐ R, for all k ∈ N and for all p ∈ N there exists a constant c 1 > 0 such that for arbitrary ε ∈ (0, 1] and for any k ∈ N 0 .
We now introduce the Colombeau algebra as the quotient For the general analysis of G(R) we refer to e.g. Oberguggenberger [Obe92].
Theorem 4.2 (Uniqueness). Assume that b is a piecewise continuous increasing function such that b ≥ b 0 for some constant b 0 > 0, and that b ′ is a positive distribu- has a unique solution u ∈ G(R + ; H s ).
Here G(R + ; H s ) stands for the space of families which are in G(R + ) with respect to t and in H s with respect to x.
Proof. Let us show that by embedding the coefficient in the corresponding Colombeau algebras the Cauchy problem has a unique solution u ∈ G(R + ; H s ). The existence follows along the lines of the proof of Theorem 2.8. Assume now that the Cauchy problem has another solution v ∈ G(R + ; H s ). At the level of representatives this means , ) ε are C ∞ (R + )-negligible being from the same equivalence class in the Colombeau algebra.
Hence, using Lemma 2.7 globally in T , we get the statement of the theorem since the function (f ε ) ε is C ∞ (R + ; H s )-negligible.
Proof of Theorem 2.9. Now we compare the classical solution u provided by Theorems 2.3 and 2.4 with the very weak solution u given by Theorem 2.9. The classical solution satisfies (4.1) For the very weak solution u, there is a representative (u ε ) ε of u such that (4.2) for a suitable embedding of the coefficient b.
. Let us rewrite (4.1) as and n ε ∈ C([0, T ]; H s−1 ). Also n ε → 0 as ε → 0 in C([0, T ]; H s−1 ). From (4.2) and (4.3) we get that ( u − u ε ) solves the Cauchy problem Since by Lemma 2.7 we have Remark 4.3. The convergence can be made on the interval [0, ∞) compared to the statement of Theorem 2.9, depending on further properties of b(t). We note that in the proof of Theorem 2.9 we apply Lemma 2.7 with To ensure the convergence of the integral in (2.13), for example under conditions of Theorem 2.4, using the estimate is finite under the conditions of Theorem 2.4. Similarly, under conditions of Theorem 2.3, using that is finite depending on further properties of b(t).

Distributional initial data case
In this section we briefly discuss the case when the Cauchy data are distributional. More precisely, we consider the following initial problem as in (2.1), but here the Cauchy data (u 0 , u 1 ) are allowed to be from H s × H s−1 with a negative s. Indeed, the main argument in allowing s to be also negative is the fact that the results of Theorems 2.3 and 2.4 could be extended as follows: where D x −σ is a Fourier multiplier with the symbol (1 + |ξ| 2 ) − σ 2 , (σ > 0); • Secondly, we consider the Cauchy problem (5.1) with v, that is • Finally, for σ > 0 we recall the Sobolev spaces Since the coefficients of the equation (5.2) do not depend on x the statements of Theorems 2.3 and 2.4 hold for any s. Thus, arguing in the same way as in the previous sections, we get the following results on very weak solutions: Theorem 5.1 (Existence). Assume that the coefficient b of the Cauchy problem (5.1) is a positive, piecewise continuous and increasing function such that b ≥ b 0 for some constant b 0 > 0, and that b ′ is a positive distribution such that b ′ ≥ b ′ 0 for some constant b ′ 0 > 0. Let s ∈ R and let the Cauchy data (u 0 , u 1 ) be in H s × H s−1 . Then the Cauchy problem (5.1) has a very weak solution of order s.
Theorem 5.2 (Consistency). Assume that b ∈ C 1 ([0, ∞)) is an increasing function such that b ≥ b 0 > 0, and that b ′ is a positive distribution such that b ′ ≥ b ′ 0 for some constant b ′ 0 > 0. Let s ∈ R, and consider the Cauchy problem (5.1) with (u 0 , u 1 ) ∈ H s ×H s−1 . Let u be a very weak solution of (5.1). Then for any regularising family b ε in Definition 2.
Then there exists an embedding of the coefficients b, b ′ into G(R + ), such that the Cauchy problem (5.1) has a unique solution u ∈ G(R + ; H s ).
We note that certain further generalisations of the obtained results are possible replacing the Laplacian with more general operators with a control on its spectral behaviour, using methods developed in [RT16a,RT17], or to dissipative wave equations on manifolds ( [DRT16]). Such questions will be addressed elsewhere.

Numerical experiments
To illustrate the theoretical results presented above, in this section we compute a sequence of solutions (u ǫ ) ǫ>0 of the regularised problem (2.4), obtained for the case of a regularisation of the synthetic piecewise continuous increasing function b(t) = 1, t < 5, 1 10 t + 3 2 , t ≥ 5. We consider the regularisation b ǫ (t) = (b * ψ ǫ )(t), t ≥ 0, of the coefficient b(t) by the convolution with the mollifier where Here C = 0.443994 so that ψ = 1. We take the initial conditions for the regularised problem (2.4) to be where δ = 0.3 and a = 0. In particular, in Figure 1, we illustrate the regularisation b ǫ (t) of the coefficient b(t) for ǫ = 0.5. To approximate the solutions of the regularised problem (2.4), we solve the equivalent system subject to the initial conditions Observe that if (p, u) is a solution of problem (6.2), then u ǫ = u is a solution of problem (2.4).
To solve numerically the initial value problem (6.2), we approximate the first derivatives p x , u x with respect to the variable x with a fourth-order finite difference scheme, and the classical fourth-order Runge-Kutta method is used for time stepping. The results for u ǫ (t, x) at time t = 60 are presented in Figure 2, for several values of the parameter ǫ. The numerical parameters used in these computer simulations were ∆x = 0.0171, ∆t = 0.0067, and the spatial computational domain was the interval [−50, 70]. In all computer simulations, we used Matlab R2016b.
In Figure 3, we display the solution u ǫ (t, x) for ǫ = 0.01 at different values of time t = 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2 before and after the interaction with the discontinuity at t = 5. Observe that the profile starts splitting into two approximately at t = 5.6, generating a second traveling wave with smaller size moving to the left.  Figure 1. In this plot, we display the function b(t) together with its regularisation b ǫ (t) obtained by convolution with the mollifier ψ ǫ (t) defined in (6.1). The regularised coefficient b ǫ was computed numerically with Matlab R2016b.
In Figure 4 we illustrate the decay with respect to time t of the solution u ǫ of the regularised problem (2.4) for some values of the parameter ǫ.
Finally, in order to illustrate the evolution of a near-singular initial pulse in problem (2.4), in Figures 5, 6, 7, we display the solution u ǫ of the regularised problem (2.4) at time t = 8, with ǫ = 0.01 and initial data in the form where f (x) = 1 π e x 2 + e 2 , with e = 0.05, e = 0.03 and e = 0.01, respectively. The coefficient b(t) and its regularisation b ǫ (t) were the same as in the previous experiments. The numerical parameters were ∆t = 0.0011, ∆x = 0.025 for the experiment with e = 0.05, and ∆t = 8E − 4, ∆x = 0.025, for e = 0.03 and ∆t = 2.28E − 4, ∆x = 0.008 for the case of the experiment with e = 0.01. In all experiments the spatial computational domain was the interval [−20, 20]. Thus, we are considering a sequence of incoming pulses that approaches to the delta function. From these pictures, we can see that the left-going traveling wave (already observed in the experiments in Figure 3) does not appear to be regularised, and instead resembles a delta function but just smaller in amplitude. We point out that this phenomenon is different from conical refraction when the splitting singularity is more regular than the original one. ǫ=0.3 ǫ=0.5 ǫ=0.7 ǫ=0.08 ǫ=0.01 6.1. Conclusions for the numerical part. The numerical experiments demonstrate that the approximation techniques work well also in the situation when the strict mathematical formulation of the problem is difficult within the classical theory of distributions. The notion of very weak solutions eliminates this difficulty yielding the well-posedness results for equations with singular coefficients. Within this approach (of very weak solutions) one can recover the expected physical properties of the equation, for example the propagation profile and the decay of the sup-norm of solutions for large times. Moreover, we seem to discover a new interesting phenomenon, presented in Figure 3 and analysed further in Figures 5-7: the appearance of a new (reflective) wave (shortly) after the singular time, travelling in the direction opposite to the main one. This can be explained as an echo effect in the original acoustic problem produced at the interfaces of discontinuity of the medium. Moreover, the reflected wave appears to be of the same regularity as the original one, just smaller in amplitude. Therefore, this phenomenon is different from the one appearing in conical refraction in the presence of multiple characteristics.