How does rescaling work with integrals?
The rescaling is nothing else than a change of variables and thus follows exactly the same rules as a change of variables.
Also, it can be conceptually helpful to recall that integration is essentially the same as the computation of an area; how do areas change under rescaling?
What numerical scheme should I use for the solution of the glacier equation?
If you are modelling a retreating glacier (more precisely, a situation where no shocks are expected to appear in the solution), you can use a finite difference scheme without any major difficulties. You need to be careful, though, with the choice of the spatial discretisation: Here it is recommended to use "upwind schemes", where the choice of the discretisation is based on the direction of the flow in the equation, in order to minimise issues with the stability of the scheme.
If your are modelling an advancing glacier (where you can expect a shock solution), a better idea is the usage of "finite volume schemes." Here, the simplest possibility is probably "Godunov's scheme."
The numerical solutions that I obtain are unstable.
If you have not done so already, change your implementation to an upwind based or finite volume based scheme. If this is what you are already using and you still observe instabilities, you might want to check the ratio between the temporal and spatial discretisation levels. Recall that all explicit numerical schemes will become unstable unless a "Courant–Friedrichs–Lewy condition" (CFL condition) is satisfied.