This is something I probably heard many times at many courses I was attending during my long years of studies. Nevertheless, I don't think I really understood what was going on physically, although I have never had problems understanding mathematics involved, i.e. equations were clear as well as transformations required to linearize them, for example to get a wave dispersion relation $\omega=\omega(k)$.
Currently, I am going (slowly) through the deformations section ("Deformations and Strains. Wave speed") in my favourite physics problem book (Savchenko et al 1981) and one problem is asking to derive an expression for the phase speed of gravity waves in a shallow water layer, assuming the wave lengths to be much larger than the depth of the water layer. Another implied condition is that the density of the fluid is constant so that we indeed focus only on gravity waves, because a fluid with varying density would allow for elastic waves (I hope this is the right term here...).
I was trying to think of perturbations and how they would be propagating and about trajectories of water parcels subjected to the above constraints (mostly mass conservation). But honestly I could not understand right away how to go about solving the problem without formally writing equations of motion and actually got some help and direction from "Essentials of Atmospheric and Oceanic Dynamics" by Geoffrey Vallis. When I finally clearly realized the path for this formal solution I wrote it down as follows.
Figure 1 shows the water layer and a perturbation propagating to the right at speed $c$.
Considering the problem to be one-dimensional with the $x$-axis pointing towards the direction of the propagation of waves, we can write the following linearized momentum and mass continuity equations.
$$ \left\{\begin{array}{l} \partial_t u &= -\frac{1}{\rho_0}\partial_x P = -g \partial_x \eta \\ \partial_t \eta &= -h \partial_x u \end{array}\right.\tag{1} $$These $u$ ($x$-component of the current) and $\eta$ (elevation relative to the unperturbed depth $h$) are small perturbations around the base state of the shallow water layer at rest. The impact of Earth rotation is ignored here, you could checkout "Essentials of Atmospheric and Oceanic Dynamics" for more details about the Earth rotation effects (i.e. Coriolis force).
Then by applying $\partial_x$ and $\partial_t$ to the first and the second equations of the system (1) respectively we obtain the following equation for $\eta$:
$$ -\frac{1}{h} \partial_{t}^2 \eta = -g\partial_{x}^2\eta \Rightarrow \\ \boxed{\frac{1}{gh} \partial_{t}^2 \eta = \partial_{x}^2\eta} $$Comparing the above to the canonical wave equation
$$ \frac{1}{c^2}\partial_t^2 \xi = \partial_{x}^2\xi + \partial_{y}^2\xi + \partial_{z}^2\xi = \Delta \xi $$we can conclude that the wave speed in our case is $c=\sqrt{gh}$. Note how it does not depend on $k=2\pi/\lambda$, i.e. the same for all wavelengths, so that the initial perturbation won't be deforming and would maintain its shape with time. This is not the case for waves in deep water.
That would close the problem for me and I would move on, but somehow I decided to dig a bit around and learned a bit more about these waves in Feynman's physics lectures, although he was considering short waves in a very deep layer (even considering capillary effects). There he describes water parcel trajectories which are not that simple as one can think. In Feynman's book and at the back of the problem book itself I found hints to the following simpler solution of this problem.
This alternative approach is visualised in Figure 2.
Assuming that we have a limit between unperturbed and perturbed parts of the fluid propagating with speed $c$ (the one we are interested to compute). Let's denote the height difference between these parts of the fluid as $\Delta h$. If we accept that water from the perturbed region would pile onto the water from the unperturbed region, we could express the water flux as follows:
$$ u(h + \Delta h) = c \Delta h $$We can simplify this expression using the assumption that $\Delta h \ll h$ and therefore:
$$ uh = c \Delta h \Rightarrow u = c \frac{\Delta h}{h} \tag{2} $$Now we can use momentum equation to get more information about propagation of the perturbed region. During time $\Delta t$ a part of the unperturbed region will gain the same speed as in the perturbed region and we can write the momentum change as follows:
$$ \Delta p = \rho_0 c \Delta t h \Delta y \cdot u \Rightarrow \frac{\Delta p}{\Delta t} = \rho_0 c h \Delta y \cdot u $$This momentum change is caused by the pressure of the perturbed region on the unperturbed region:
$$ F = \rho_0 g \Delta h \cdot \Delta y h = \rho_0 c h \Delta y \cdot u $$Then simplifying we get
$$ g \Delta h = c \cdot u $$And substituting $u$ from equation (2) we get the following:
$$ g \Delta h = c \cdot c \frac{\Delta h}{h} \Rightarrow c=\sqrt{gh} $$