Why is normal ordering even a valid operation in the first place? I mean it can give us some nice results, but why can we do the ordering for the operators like that?
Is its definition motivated by the relation between the normal ordering and the time-ordered product, which is basically the content of Wick's theorem?
In classical physics, quantities are ordinary, commuting $c$-numbers. The order in which we write terms in expressions is of no consequence. In quantum field theory (QFT), on the other hand, quantities are described by operators that, in general, don't commute.
Classical physics is a low-energy approximation of quantum physics - the road from quantum to classical physics ought to be unambiguous - and this is way the way nature goes, from high to low energies. The inverse - the road from classical to quantum, that we take to try and reconstruct the high-energy physics - however, is ambiguous, because of ordering ambiguities in non-commuting quantities.
When we normal order expressions after canonical quantization, we are correcting those ambiguities.
This occurs for the zero-point energy in the Hamiltonian $$ H = \int \frac{d^3p}{(2\pi)^3} E_p \left(a_p^\dagger a_p + \frac12[a_p^\dagger,a_p]\right) $$ You might hear it argued that the since the vacuum energy is unobservable, we are free to throw away the divergent piece (the commutator). Such an argument doesn't work for the charge operator, $$ Q = \int \frac{d^3p}{(2\pi)^3} E_p \left(a_p^\dagger a_p - b_p b_p^\dagger\right) $$ A charged vacuum would have observable effects. The best argument for normal ordering is that it is a rule for removing ordering ambiguities that results in e.g. a neutral vacuum.
Ordering ambiguties also occur in general relativity, when one promotes commuting ordinary derivatives $\partial_\mu$ to non-commuting covariant derivatives $\nabla_\mu$.
The choice of normal ordering prescription $:~:$ is usually adjusted to the choices of bra vacuum state $\langle \Omega|$ and ket vacuum state $|\Omega\rangle$, so that
$$\langle \Omega|:\hat{\cal O}_1\ldots \hat{\cal O}_{n>0} : |\Omega\rangle~~=~~0 , \qquad\qquad \langle \Omega|\Omega\rangle~~=~~1 .$$
The relation of normal ordering prescription to Wick theorem and other operator ordering prescriptions are typically only secondary. For motivation, see also e.g. this Phys.SE post.
normal ordering is a valid operation provided one can undo it by an appropriate choice of counterterms (of existing couplings or field renormalisations). (How this is done in practice is explained here: http://arxiv.org/abs/1512.02604.)
As other answers mention, it was originally (in QED) about getting a neutral vacuum. It is useful to go back to Schwinger's old version of QED, before Dyson's approach became accepted. See Pauli: Selected topics in field quantization.
Pauli presents both ways of looking at it: 1) define the electric current as sum of two terms (p.20 [6.4]), such that the vacuum expectation number operators give zero [6.7-6.9]; and later see that 2) normal ordering, i.e. dropping the infinite vacuum constant, has the same effect when taking the ordinary definition of the current, as far as the S-matrix is concerned (see p.141: "Not grouping together terms with the same argument is a substitute for the neglected subtraction of the vacuum current...").
This is how the 'negative energy sea' was eliminated. There is nothing incorrect about it, other than that the true content of QFT, i.e. renormalization, reaffirms the S-matrix point of view, 2).
The first approach tried to work with operators without reference to the S-matrix, as long as possible, whereas the second concedes from the start that only S-matrix results will be useful.
(Caution: Pauli's book is interesting but not easy to read.)
Here is a simple non-physics explanation of why the operation is useful:
Consider the equation:
$$f(x,y)e^{xy} = :f(x,\partial_x): e^{xy}$$
The equation only works if the operator function $f$ on the right is 'normal-ordered' with all the derivatives on the right of the terms.
So if we ask the question, what operator corresponds to the function $f$ that makes this equation work? The answer is the 'normal-ordered' operator made from $f$. So we can see that `normal ordering' does make have a practical use.
