16 months after asking the question, I've run into a different and quite physicsy answer, which I hope is useful to someone.
Suppose that in order to determine $m$ parameters of some model:
$$\text{Output} = f(\text{Inputs, parameters})$$
we have conducted $N>m$ experiments. We want to use the information from these experiments to best choose the $m$ parameters (so that the output of the model and the actual experimental value are as close together as possible).
Now comes the physicsy part: lets construct an N-dimensional phase space so that the $N$ (experimentally determined) outputs from our $N$ experiments are represented by a single point in this space (the Cartesian coordinates of this point are the outputs from each experiment). Call this the 'data-point'.
Secondly, if we choose an arbitrary set of parameters for our model, we can use the inputs for each experiment to construct a 'predicted output' for each experiment (by parsing the inputs through our model). There will be $N$ predicted outputs (one for each experiment) and these form a second point in the phase space, say the 'prediction-point'. As we varying the parameters this point moves about in an $m$-dimensional subspace of the phase space. And this is the important point:
The sum of the squares of the error terms (SSE) is the square of the distance between these two points in the $N$-dimensional phase space, just by Pythagoras' Theorem.
So minimising the sum of squares error is equivalent to minimising the distance between the data-point and the prediction-point in the $N$-dimensional phase space - a very natural way of calibrating our model.
Finally, from this Gauss' result makes some sense - if the data point is allowed to vary normally with 0 mean and equal variances, the error will be spherically-symmetric around the data-point, and so the closer our prediction-point is to the data-point in space, the better, and minimising this distance should give the maximum-likelihood estimator.