Source terms, such as those due to bathymetry in the shallow water equations, need to be integrated in a special way in order to preserve physical steady states. Is there a general way to construct well-balanced methods, or does it require special techniques for each equation?
The short answer is: it requires specific work for different equations, but there are some general techniques that suggest how to do it. Essentially, given a first order evolution PDE
$$u_t = Au + Bu$$
where $A,B$ are some (possibly differential) operators, the steady states are those for which
$$Au + Bu=0.$$
It is common to use a splitting approach in which $A$ and $B$ are discretized in different ways. Then there will be truncation errors associated with each of these discretizations, and the truncation errors will generally not cancel even in the case of a steady state. A classical example (as mentioned in the question) is the shallow water equations with bathymetry, in which $A$ represents convective terms and $B$ represents the momentum forcing due to the variable bottom height. There are many papers published in the last few years that give different ways to exactly maintain the steady state solutions.
One approach that I like is the use of f-wave Riemann solvers as proposed by Bale et. al.. The idea is to discretize the convective terms with a Godunov-type method, but to subtract off the contribution from the other terms inside the Riemann solver. Then in the case of the steady state, no waves are generated. However, this requires that the convective and source terms be calculated exactly (so as to cancel exactly). That's possible for the shallow water equations, but more difficult for many other systems.