Modflow comes to mind as a modelling tool. Note that a lot of things influence this process. Initial state, to what extent is the soil already saturated, where is the groundwater table? What is the hydraulic conductivity of the subsurface? Vegetation, surface temperature, wind, humidity, incoming radiation all influence evapotranspiration losses too.

Approach 1: create a density model of the planet and integrate(sum) this model over all volume elements assuming every elements contributes as a point mass. For the Earth, the largest non-spherical anomaly comes from its flattened shape which cause satellite to deviate from a pure Kepler orbit. The further away, the more the planet looks like a point mass, and the more satellite orbits approach a Kepler orbit (perfect ellipse).

Approach 2: use satellite orbit or astronomy observations in combination with our knowledge that the Earth's (or other planet) external gravity field obeys the Laplace equation. This allows a gravitational model to be estimated for the external field. However, knowing the external field is not enough to know the internal mass distribution (infinite possibilities still generate the same external gravity field)

Additional contributions from other far objects (moons, other planets) can generally added as additional point masses, and are considered tidal accelerations. This requires knowledge of the position of those objects relative to the planet, which can be taken from a so-called planetary ephemeris.