You don't seem to be considering interactions here, so you don't even need the GPE. This is just about Bose-Einstein condensation (non-interacting phenomenon).
Without interactions, number operators are defined and the state is a Fock state:
$$ |\Psi\rangle = \left ( \prod_{j=1}^\nu \frac{1}{\sqrt{N_j!} }(a_j^\dagger)^{N_j} \right ) |0\rangle, $$
where $|0\rangle$ is the vacuum.
For BECs, you generally look at the ground state.
The expectation value of the number operator in the ground state is:
$$ \langle n \rangle = \langle \psi^\dagger\psi\rangle = \langle 0| \psi^\dagger \psi |0\rangle \propto \sum_k \langle 0|a^\dagger_k a_k|0\rangle = 0, $$
because $a_k |0\rangle =0$. Just like in quantum mechanics (first quantisation). So the number operator is quite useless for the ground state, i.e. for a BEC.
Now.
You know that a BEC breaks the $U(1)$ symmetry corresponding to particle conservation. Hence,
for $T>T_{\mathrm{c}}$, the number operator makes sense and you can build a Fock state like above. Hence, $\langle \hat\psi \rangle = \langle 0 | \hat\psi|0 \rangle \propto a_0 |0\rangle = 0.$
for $T<T_{\mathrm{c}}$, the number operator is not useful any longer as you've broken $U(1)$ and hence particle conservation -- you went from a complex field to a real field, $\hat\psi \in \mathbb{C} \rightarrow \hat\psi \in \mathbb{R}$.
Hence:
$$\langle \hat\psi \rangle = \langle 0 | \hat\psi|0 \rangle \neq 0,$$
which is one of the definitions for the BEC phase.
In fact, the BEC is actually a coherent state of the $|\mathbf{k}=0\rangle$ mode:
$$ |BEC\rangle = \mathrm{e}^{-N/2} \mathrm{e}^{\sqrt{N}a_0^\dagger} |0\rangle, $$
where $|0\rangle$ is, again, the vacuum, and $N$ is the expectation value of the number operator, since the state is not a Fock state and hence does not have a fixed eigenvalue for $N$. In fact, you know that a coherent state has a fixed phase $\theta$, and the coherent state $\leftrightarrow$ Fock state correspond to the extremes of the commutation relation $[\hat N, \hat \theta] \propto \mathrm{i} \hbar$ (with an associated uncertainty principle).
Ok, your ground state is now $|BEC\rangle$. So let's use this for the averages from earlier:
$\langle BEC | \hat \psi \psi | BEC \rangle = \langle n \rangle = N$ (mean number of atoms),
$\sigma_n = ... = \sqrt{N}$ (standard deviation of atom number: non-zero because it is not a Fock state,
$\langle BEC | \hat \psi(\mathbf{r},t) | BEC \rangle = \langle BEC| \sum_k \phi_k(\mathbf{r},t) \hat{a}_k |BEC\rangle = \phi_0(\mathbf{r},t)\langle BEC | BEC \rangle = \phi_0(\mathbf{r},t)$,
because a coherent state $|BEC\rangle$ is an eigenstate of the annihilation operator $a_k$, and the BEC is only built of $a_0$.
Now since $\psi^\dagger$ creates an excitation of the field $\psi$ at some position and some time, $\phi(\mathbf{r},t)$ is basically the spatial profile of the excitation.
This whole thing was done in second quantisation, so it's a bit weird to jump to first quantisation and call this "a wavefunction". It still has a "normalisation" interpretation in that the state is a 100% pure BEC, so you have a 100% probability of finding anything at $|\mathbf{k}=0\rangle$. However, the concept of "local particle density" is a bit ill-defined since this is not a Fock state and particle number $N$ has an uncertainty $\sqrt{N}$. Of course, for large number of particles $\sqrt{N}/N$ becomes negligible so you can just treat it as a "wavefunction" and not care about the nuanced meaning of particle number.
