Throughout, let's assume that the ground state energy of the system under consideration is zero.
Chay Paterson has addressed your question in the case of a gas of bosons in which the number of particles is not conserved, but from the wording of your question, it seems that you're concerned about the case in which the total number of particles is fixed.
For a system of bosons with a fixed number $N$ of particles, the answer to your question is
No. The chemical potential is nonzero for all $T>0$.
You point out, however, that it is often said that below the critical temperature $T_c$, the chemical potential is zero, so what's going on? The resolution is essentially that
the chemical potential very well-approximated by zero for almost all temperatures below the critical temperature, but it is never exactly zero.
As you cool the system down from above the critical temperature, the number $N_e$ of particles in excited states get's lower and lower. On the other hand, the number of particles in the excited states at any given temperature is bounded above, and this bound decreases as a function of temperature like $T^{3/2}$ (for non-relativistic systems in three dimensions). At a particular sufficiently low temperature (the critical temperature), this upper bound is lower than the total number of particles in the system and particles are forced into the ground state. But at this point, the chemical potential does not drop exactly to zero. It does, however, get very small very quickly as the temperature decreases and the number of particles in the ground state increases.
In fact, If we look at the number of particles in the ground state as a function of temperature
$$
N_0 = \frac{1}{e^{-\mu(T)/kT}-1}
$$
which gives
$$
\mu(T) = -kT\ln\left(1+\frac{1}{N_0}\right)
$$
then we see that the chemical potential is always strictly less than zero because the argument of the log is always strictly greater than $1$. But as you dial the temperature down below the critical temperature, and as the number of particles in the ground state increases towards $N$, the total number of particles in the system, which is presumably very large, and the chemical potential decreases since the argument of the log approaches $1$.
It feels like you are going to fast in your way to think.
The difference between the chemical potential and the ground state energy is clear :
The ground-state energy of your system is here $\epsilon_\textbf{k}=0$ correponding to the ground-state $|\textbf{k=0}\rangle$, which is macroscopically occupied in a BEC (i.e. $N\sim N_0$).
The chemical potential $\mu$ is the increment in energy that you give to your condensate by adding one particle. It is generally defined using free energy $F$ :
$$\mu=\left(\frac{\partial F}{\partial N}\right)_{V,\;T}$$
If you want to evaluate $\mu$, considere the ground-state occupation :
$$N_0=\frac{1}{e^{-\beta\mu}-1}$$
If you want $N_0$ to be the biggest possible (macroscopic occupation), you must fulfill the condition $e^{-\beta\mu}\rightarrow 1$, in other words $\mu\rightarrow 0$.
For larger resolution on $\mu$, one can do a Taylor expansion on the expression of $N_0$ :
$$N_0\underset{\mu\rightarrow 0}{\sim}-\frac{1}{\beta\mu}$$
Then,
$$\mu\simeq -\frac{1}{N_0\beta}=-\frac{k_\mathrm{B} T}{N_0}=-\frac{k_\mathrm{B} T}{N-N_\mathrm{exc}}$$
Best Answer
That "first"/"then" isn't the logic used in deriving the Bose-Einstein condensate, but it is an intuitive picture of what's going on.
So, we have for our grand canonical potential something like: $\Omega=T \sum_{i} \log(1-e^{\beta(\mu-\varepsilon_i)})$, and we can find the number of particles through $N=-\frac{\partial \Omega}{\partial \mu}$. Those are exact, true equations of quantum statistical mechanics. Here is where your confusion lies:
We approximate $\Omega=T \int \nu(\varepsilon)\log(1-e^{\beta(\mu-\varepsilon)})d\varepsilon$, where $\nu$ is the density of states. This is a good approximation in the Fermi gas, where the low energy states are filled but only by a few out of $N\approx 10^{23}$ particles. This is a good approximation for an ideal gas, where $T$ is large and most occupation is in higher energy states. But this is a terrible approximation for BECs because the behavior of the low energy states isn't captured at all by the integral!
So what are our options now? We probably don't want to work directly with the sum, but we can't just take an integral. Some books like to pull out the ground state term and approximate the rest with an integral, and I'll do that here.
$$\Omega=T \log(1-e^{\beta \mu})+\int \nu(\varepsilon)\log(1-e^{\beta(\mu-\varepsilon)})d\varepsilon$$
This gives... (denote $z=e^{\beta \mu}$)
$$N=\frac{z}{1-z}+\int \nu(\varepsilon)\frac{1}{z^{-1} e^{\beta \varepsilon}-1}d\varepsilon$$
For low temperature, the rightmost term is basically constant and the leftmost one varies dramatically ($\frac{z}{1-z}\to\infty$ as $\mu \to 0$), allowing you to match any desired value of $N$.
Now, you could ask the question: At a finite temperature $T$ with particle number $N$, if I add one more to $N$, how much is split between the ground state and how much is split between the excited states? Well, if $|\mu|$ is small (this is equivalent to $T<T_c$), $z\approx 1-\beta |\mu|$, $\nu=C \varepsilon^{1/2}$, and $\int \frac{\varepsilon^{1/2}}{z^{-1}e^{\beta\varepsilon}-1}d\varepsilon=T^{3/2} \Gamma(3/2) \mathrm{Li}_{3/2}(z)$. Mathematica tells me $\mathrm{Li}_{3/2}(1-\beta |\mu|)\approx \zeta(3/2)-2\sqrt{\pi \beta |\mu|}+O(\beta \mu)$, so overall:
\begin{align*} N&\approx\frac{z}{1-z}+\int \nu(\varepsilon)\frac{1}{z^{-1} e^{\beta \varepsilon}-1}d\varepsilon\\ &= \frac{1}{\beta |\mu|}-1+C T^{3/2} \Gamma(3/2)\left( \zeta(3/2)-2\sqrt{\pi \beta |\mu|} \right)+O(\mu) \end{align*}
For large $N$ and small $\mu$, the contribution of the first term, $\frac{1}{\mu}$, dominates by far! It tells you that "as $N\to \infty$, the ground state starts to hold almost all the states".
The contribution of the proceeding terms include a constant term, which is the notion that you fill excited states first, and a correction to the notion, which says you leave a tiny bit of the excited states unfilled, proportional to $\sqrt{|\mu|}$.
You can go on and on with this stuff at varying degrees of rigor. Pathria Statistical Mechanics appendix F does so in a way that doesn't pull the first term out of the integral, but instead works with the sums directly by writing $\varepsilon$ in terms of the energy levels of a particle in a 3D box.
[edit/footnote]: note that the integrals I take are from $\varepsilon=0$ to $\varepsilon=\infty$. I assume the ground state to have zero energy! This is fine, because all we care about is the difference $\mu-\varepsilon_i$. I can shift $\varepsilon$ by whatever energy I want, so long as I shift $\mu$ correspondingly to keep the difference the same. The choice $\varepsilon_0=0$ gives the usual $\mu\to 0$ as $T\to 0$ law, but if we shifted everything to $\varepsilon_0=10 \mathrm{ Joules}$ or whatever, you would certainly have $\mu\to 10 \mathrm{ Joules}$ as $T\to 0$. The equations would be a lot uglier but the physics would be the same.