Preliminary Note: The wave equation you present is usually used for quantum mechanics, where electromagnetic forces are significant. On the other hand, your recommended examples deal with mechanical waves. Throughout this question, I will be assuming we are dealing with fluid acoustics, since I am most familiar with that and because you seem to be requesting something a litte more macroscopic.
We have to start somewhere, and I will start assuming you know the fundamental equations of acoustics, in particular the inviscid and adiabatic equation of state and the conservation equations for mass and momentum. There are many texts that derive these equations; one example is Fundamentals of Physical Acoustics by Blackstock (2000). These equations may be linearized and written for a bulk medium, respectively, as
\begin{gather}
p = c_0^2\rho', \\
\rho_0\nabla\cdot\vec v + \frac{\partial\rho'}{\partial t} = 0, \\
\nabla p + \rho_0\frac{\partial\vec v}{\partial t} = 0,
\end{gather}
where $p$ is the acoustic pressure, $c_0$ is the small-signal sound speed, $\rho'$ is the acoustic mass density variation, $\vec v$ is the acoustic particle velocity, $\rho_0$ is the ambient density, $t$ is the time, and $\nabla$ is the gradient operator. Physically, the $\partial\rho'/\partial t$ term is the rate of increase of the mass density, and it is related to the rate at which the fluid is being compressed, the volumetric strain rate $\nabla\cdot\vec v$.
If mass were being injected into the system (think of the pebble falling into the water), then there would need to be an additional term in the conservation of mass equation. If we assume that the mass is injected at a point $\vec x'$, then we may write the mass conservation equation as
\begin{equation}
\rho_0\nabla\cdot\vec v + \frac{\partial\rho'}{\partial t} = \delta(\vec x-\vec x')\frac{\partial m}{\partial t},
\end{equation}
where $m$ is the mass being injected as a function of time. A similar argument coule be used to explain how a force being applied at the point would lead to an additional term to the momentum conservation equation, which leads to
\begin{equation}
\nabla p + \rho_0\frac{\partial\vec v}{\partial t} = \delta(\vec x-\vec x')\vec f,
\end{equation}
where $\vec f$ is the overall force being applied.
These expressions are for point sources. However, if we re-interpret $m$ and $\vec f$ as mass and force distributions that vary over space, we can drop the $\delta$ functions. Substituting in the equation of state we then obtain
\begin{gather}
\rho_0\nabla\cdot\vec v + \frac{1}{c_0^2}\frac{\partial p}{\partial t} = \frac{\partial m}{\partial t}, \\
\nabla p + \rho_0\frac{\partial\vec v}{\partial t} = \vec f.
\end{gather}
With these two equations we may eliminate $\vec v$. To do this, we take the time derivative of the first equation and the divergence of the second equation, and then subtracting the first from the second. This process leads to
\begin{equation}
\nabla^2 p - \frac{1}{c_0^2}\frac{\partial^2 p}{\partial t^2} = \nabla\cdot\vec f - \frac{\partial^2m}{\partial t^2}.
\end{equation}
Using the d'Alembertian operator ($\square=\nabla^2 - \partial^2/\partial t^2$) as you did in your post, we may then write
\begin{equation}
\square p = S,
\end{equation}
where $S$ is the source term given by
\begin{equation}
S = \nabla\cdot\vec f - \frac{\partial^2m}{\partial t^2}.
\end{equation}
For acoustics we don't typically invoke the potential term as you wrote it; that form is usually used in quantum mechanics to account for the electromagnetic (or other) force terms. A potential-like term does appear if you account for gravity in the sound systems, but except in cases such as infrasound (very long wavelengths) the gravitational terms do not contribute significantly. Without deriving it [the derivation may be found in Acoustics by Pierce (1989)], I will give the result from assuming an equilibrium atmosphere with gravity and only allowing motion in the vertical, $z$, direction:
$$ \frac{\partial^2\tilde p}{\partial z^2} - \frac{1}{c_0^2}\frac{\partial^2\tilde p}{\partial t^2} - \frac{\gamma g}{2c_0^3}\tilde p = 0, $$
where $\gamma$ is the ratio of specific heats and $\tilde p=p/\sqrt{\rho_0}$. Thus, using the d'Alembertial operator again (but only in the $z$ direction), we may write
$$ \square\tilde p + V\tilde p = 0, $$
where the potential term is
$$ V = -\frac{\gamma g}{2c_0^3}. $$
Including both gravity and a mass term would lead to an equation that uses both the potential and source terms. It also can get a lot more complicated with many additional terms as more interesting situations and physical phenomena are accounted for. The world is fascinating, and sometimes simple equations are not sufficient to describe it all.
