From what I have gathered, the Fermi level of a semiconductor is equal to the chemical potential of the electrons, i.e., it's the work required to add an electron to the system.
We've got to be a bit careful here. The chemical potential is the change in internal energy of the system when we add one electron at constant entropy. More generally, we would have that
$$\mu = \Delta U - T\Delta S$$
Only at $T=0$ do we have that $\mu = \Delta U$ when a single particle is added. This suggests that $\mu(T\rightarrow 0) = E_c$, the energy at the bottom of the conduction band, which turns out to be true in contrast to the conventional wisdom which holds that at $T=0$, the Fermi level lies in the middle of the gap.
To understand how this misconception arose and why it persists, we need to re-examine the derivation of the Fermi-Dirac distribution function
$$f(E) = \frac{1}{e^{(E-\mu)/kT}+1}$$
which gives the probability that an energy level with energy $E$ is occupied. I will follow the derivation given in ref 1.
Adopting the canonical ensemble perspective, we fix the temperature $T$ and particle number $N$; a microstate then consists of a list $\{n\}\equiv\{n_1,n_2,\ldots\}$ of the occupation numbers of each single-particle energy level.
The partition function becomes
$$Z(T,N) = \sum_{\{n\}\in M_N} e^{-\beta E\big(\{n\}\big)} = \sum_{\{n\}} \exp\left[-\beta\sum_i n_i \epsilon_i\right]=\sum_{\{n\}\in M_N} \prod_i e^{-\beta n_i \epsilon_i}$$
where $M_N$ is the set of microstates such that $\sum_i n_i = N$, the summation over $i$ corresponds to the single-particle states, and $\epsilon_i$ is the energy of the $i^{th}$ state. This summation cannot be easily performed for arbitrary $N$ because of the constraint $\{n\}\in M_N$. However, we can extract the distribution function as follows.
Consider in particular the $I^{th}$ state. The probability that the state is occupied is given by
$$f_I(T,N) = \sum_{\matrix{\{n\}\in M_N\\n_I=1}}e^{-\beta E\big(\{n\}\big)}/Z(T,N)$$
where the sum in the numerator is taken over all microstates $\{n\}\in M_N$ in which $n_I = 1$. Since either a state is occupied or it is not, we may equivalently write this as $1$ minus the probability that the state $I$ is not occupied, i.e.
$$f_I(T,N) = 1-\sum_{\matrix{\{n\}\in M_N\\n_I = 0}} e^{-\beta E\big(\{n\}\big)}/Z(T,N)$$
Next, we note that each $N$-particle state with energy $E$ in which state $I$ is not occupied corresponds to a unique $N+1$-particle state with energy $E+\epsilon_I$ in which state $I$ is occupied and vice-versa. Therefore, we may write
$$f_I(T,N) = 1- \sum_{\matrix{\{n\}\in M_{N+1} \\ n_I = 1}} e^{-\beta\bigg(E\big(\{n\}\big)-\epsilon_I\bigg)} / Z(T,N)= 1-e^{\beta \epsilon_I} f_I(T,N+1) \frac{Z(T,N+1)}{Z(T,N)}$$
Recalling that $Z(T,N) = e^{-\beta \mathcal F(T,N)}$ with $\mathcal F$ the Helmholtz potential, we observe that $$Z(T,N+1)/Z(T,N)=\exp\bigg(-\beta\big(\mathcal F(T,N+1)-\mathcal F(T,N)\big)\bigg)\equiv \exp\big(-\beta \mu(T,N)\big)$$ and finally obtain
$$f_I(T,N) = 1 - e^{\beta(\epsilon_I-\mu)}f_I(T,N+1)$$
where the chemical potential has been defined as $\mu(T,N)\equiv \mathcal F(T,N+1)-\mathcal F(T,N)$.
The Fermi-Dirac distribution is easily obtained under the assumption that $f_I(T,N)\approx f_I(T,N+1)$. The conventional wisdom is that this assumption is valid universally, because it is seemingly ridiculous to assume that the addition of another electron to a system cannot have any meaningful effect on the occupation probabilities. However, this is not true for a band insulator at $T=0$.
Let $N_0$ be the number of electrons which would completely fill the valance band and leave the conduction band empty, and let $I$ be the lowest energy state in the conduction band. Clearly $f_I(0,N_0)= 0$, but the addition of a single additional electron yields $f_I(0,N_0+1)=1$, rendering our assumption invalid.
A correct treatment of this problem can be found in ref 2. The result is that the low-temperature behavior of $\mu$ is given by
$$\mu= E_c -\frac{1}{2}\Delta+ kT \ln\left[\left(\frac{m_v}{m_c}\right)^{3/4} \alpha(T,V)\right] -kT N_i(T,V) \left[\alpha(T,V)+\alpha^{-1}(T,V)-2\right]$$
$$N_i(T,V) = \frac{1}{4}V \left(\frac{2kT}{\pi\hbar^2}\right)^{3/2} (m_vm_c)^{3/4} e^{-\Delta/2kT}, \qquad \alpha(T,V) = \frac{1+\sqrt{1+4N_i^2(T,V)}}{2N_i(T,V)}$$
where $E_c$ is the energy at the bottom of the conduction band, $\Delta$ is the band gap, and $m_v$ and $m_c$ are the effective masses at the top of the valence band and bottom of the conduction band, respectively. This rather unpleasant expression yields the following behavior as $T\rightarrow 0$ (figure taken from ref 2):

Ashcroft and Mermin, p.40-42
M. R. A. Shegelski, Solid State Commun. 58, 351–354, 1986