Here is a very nice animation showing the motion of water molecules under the effect of a water wave, which is taken from this site. For a hardcore mathematical treatment of this fascination subject check out this book.
Regardless of the above complications it is however a matter of observation that the height of the water surface indeed behaves like a transverse wave, and I shall therefore assume that it satisfies the two-dimensional wave equation:
\begin{equation}
\frac{\partial^2 h }{\partial x^2}
+
\frac{\partial^2 h }{\partial y^2}=
\frac{1}{v^2}\frac{\partial^2 h }{\partial t^2}
\end{equation}
I am interested in examining cylindrically symmetric solutions of the above equation. Before doing so it is useful to take a look at solutions of the one-dimensional wave equation and spherically symmetric solutions of the three dimensional wave equation. The solutions for these cases have a very nice and simple form and are intuitively simple to understand. So, in one dimension, for any arbitrary functions $f$ and $g$, $u(x,t) = f(x–vt)+g(x+vt)$ is a solution of the wave equation. This can be easily verified by substituting it in the one dimensional wave equation.
For spherically symmetric solutions in the three-dimensional case, the wave equation reads
\begin{equation}
\frac{1}{r^2}\frac{\partial }{\partial r} \left( r^2 \frac{\partial h}{\partial r} \right)=
\frac{1}{v^2}\frac{\partial^2 h }{\partial t^2}
\end{equation}
By direct substitution it can be verified that $(1/r)[f(r-vt) + g(r+vt)]$ is the general solution of the above equation, in other words, the solutions are spherically symmetric waves traveling either towards the center or away from the center, with the amplitude proportional to $1/r$.
Hence in one dimension, the amplitude is independent of $x$, and in three dimensions the amplitude varies as $1/r$. Let us understand this physically. For simplicity let us assume that the functions $f$ and $g$ are sinusoidally varying functions. Let us further assume that the medium in which the wave travels consists of an infinite number of simple harmonic oscillators oscillating as the wave travels. In one dimension let the number of such oscillators per unit length be a constant $\lambda$. A small section of the rope of length $dx$ will thus have $\lambda dx$ oscillators. The energy of a simple harmonic oscillator is proportional to the square of its amplitude. Let the amplitude of each oscillator in this section $dx$ be $A(x)$. Since there are $\lambda dx$ oscillators, the total energy in length $dx$ is proportional to $λdxA(x)^2$. This should be independent of $x$ (This can be understood in various ways: either (1) by seeing that a sinusoidal function is linearly homogeneous, or (2) by physically assuming that someone is physically moving one end of the rope thus pumping in energy at a constant rate, so that in a section $dx$ of the rope, the energy that enters from the left must be equal to the energy that leaves from the right, otherwise there will be sinks or sources of energy, or (3) by considering a single pulse traveling across the rope and seeing that the number of oscillators per unit length is independent of $x$.) Since $\lambda$ is independent of $x$, it follows that $A(x)$ also must be independent of $x$.
What happens in the spherically symmetric three dimensional case? Well, let us assume that the number of oscillators per unit volume in space is a constant $\mu$. Consider a thin spherical shell of thickness $dr$ a distance $r$ away from the origin. The volume of this shell is $4\pi r^2dr$, hence the number of oscillators in this shell is $4\pi\mu r^2dr$. Let the amplitude of each oscillator be $A(r)$, implying the energy of each oscillator is proportional to $A(r)^2$. Multiplying this with the total number of oscillators in this shell gives us the total energy in this shell, which is thus proportional to $A(r)^2\cdot 4\pi\mu r^2dr$. Again, this must be independent of $r$, which is possible only if $A(r)$ is proportional to $1/r$.
Using a similar logic as above, one could conclude that the general solution for the cylindrically symmetric wave equation in two dimensions is $\rho^{-1/2} f(\rho–vt)$, where $\rho$ is the distance of a point from the $z$-axis in cylindrical coordinates. Indeed, a thin circular disc of thickness $d\rho$ and a distance $\rho$ away from the origin will have $\alpha\cdot 2\pi\rho d\rho$ number of oscillators, where $\alpha$ is the number of oscillators per unit area. If $A(\rho)$ is the amplitude of each oscillator in this disc, the total energy is proportional to $A(\rho)^2α2π\rho d\rho$, which again must be independent of $\rho$, suggesting that $A(\rho)$ varies as $\rho^{-1/2}$.
Unfortunately, this answer is wrong. This can be seen by explicitly substituting $\rho^{-1/2} f(\rho–vt)$ in the cylindrically symmetric two dimensional wave equation:
\begin{equation}
\frac{\partial^2 h }{\partial \rho^2}
+
\frac{1}{\rho}\frac{\partial h }{\partial \rho}=
\frac{1}{v^2}\frac{\partial^2 h }{\partial t^2}
\end{equation}
This equation is not satisfied for the above form of the wave function. Indeed, there exists no $n$ for which $ρ^n f(ρ–vt)$ will satisfy the above equation. For, if you substitute
\begin{equation}
h(\rho,t)=\rho^n f(\rho-vt)
\end{equation}
in the above cylindrically symmetric two dimensional wave equation, an explicit calculation for the left hand side gives:
\begin{equation}
\frac{\partial^2 h }{\partial \rho^2}
+
\frac{1}{\rho}\frac{\partial h }{\partial \rho}
= n^2\rho^{n-2}f + (2n+1)\rho^{n-1}f' + \rho^n f''
\end{equation}
where $f'=\left. \frac{df(x)}{dx}\right|_{x=\rho-vt}$, $f''=\left. \frac{d^2f(x)}{dx^2}\right|_{x=\rho-vt}$, etc. This can satisfy the wave equation only when the first two terms (i.e. the $f$ and $f’$ terms above) vanish. The second term can indeed be removed by choosing $n = -1/2$, but the first term remains. (In contrast, for three dimensions the choice $n = -1$ satisfies the wave equation by eliminating the first two terms in the corresponding equation.)
However, for large $\rho$, the above guess does satisfy the wave equation, since in the large $\rho$ limit the first term also goes to zero. Hence the guess $\rho^{-1/2} f(\rho–vt)$ is strictly speaking a solution only in the large $\rho$ limit. The correct solution is given by the so-called Bessel functions, which indeed show an asymptotic $\rho^{-1/2}$ behavior in the large $\rho$ limit. (As an exercise one could try out the ansatz $g(\rho)f(\rho-vt)$ to see if an exact form for the amplitude $g$ can be found. This is physically more intuitive than the separation of variables method usually used.)
In our case we can always choose the boundaries at which reflection occurs to be so far away from the origin that the above asymptotic form for the amplitude holds.
Understanding multiple reflections of water waves in a square boundary¶
Having gained a bit of clarity regarding reflection of wave at a single boundary and the nature of the solution of the wave equation in one, two and three dimensions, it is time to come back to our original aim, viz. to understand multiple reflection of waves within an enclosed two dimensional region.
For simplicity, I decided to focus on a water surface enclosed in a square region and trace the multiple reflections of a single circular pulse starting from the center. From our discussion, the pulse, which is described by $f(\rho–vt)$, travels out at the speed $v$, and its amplitude falls off as $\rho^{-1/2}$. The function $\rho^{-1/2}f(\rho-vt)$ describes the height of the water above equilibrium at distance $\rho$ from the origin and at time $t$. I chose a Lorentzian form for $f(\rho-vt)$, although a Gaussian would have done as well. Thus if there are no boundaries, the height of the water above equilibrium at distance $\rho$ from origin and time $t$ is:
\begin{equation}
h(\rho,t) = \frac{1}{\sqrt{\rho}} \cdot \frac{\eta}{(\rho-vt)^2+\eta^2}
\end{equation}
where $\eta$ fixes the width of the pulse. The above just describes a circular pulse moving with velocity $v$ away from the origin. In the following animation the pulse is shown as a circle moving away from the origin.