Fermi Level: chemical potential vs. electrochemical potential In solid-state physics it is understood that the Fermi-Level is the electrochemical potential. 
The Fermi-Level is defined to:

However, in thermodynamics this formula is referred to as the chemical potential. [e.g. G Job and F Herrmann, Chemical potential—a quantity in
search of recognition, Eur. J. Phys. 27 (2006) 353–371, Equation (10)]
In literature on solid state physics [e.g. Sze], I can not find a reasoning why the Fermi-Level should also reflect the electric potential in addition to the chemical potential.
For the electric potential I would expect a solution to Poisson's Equation to be part of the definition for the Fermi-Level, which obviously is not the case. 
I am wondering if there is proof for the Fermi-Level is the electrochemical potential, and not only the chemical potential. 
 A: In the Fermi-Dirac distribution, the ``$\mu$'' is actually the
total chemical potential (ie the electrochemical potential in our
case) - but one should also take into account the influence of the electrical potential on the density of states. To see how and why, it is interesting to compare two situations, with and without electrical potential.
Consider first a reference situation, where the energy of particles
depends on the wave-vector ${\bf k}$ and the band index $n$ through
the dispersion relation $\epsilon_{n,{\bf k}}$. The corresponding
density of state is $D_{0}(\epsilon)$. Focusing on the conduction
band for instance, the density of state is 0 for energies below $E_{C}$,
and evolves as a square root above this threshold. The chemical is
the same as the electrochemical potential, and is uniform through
the system. The density of particles is also uniform.
Consider now that an electrical potential $\phi({\bf r})$ is added
to the system, such that particles now have a total energy $\epsilon_{{\rm tot}}=\epsilon_{n,{\bf k}}-e\phi({\bf r})$.
The density of state is shifted as compared to the previous situation
$D(\epsilon_{{\rm tot}},{\bf r})=D_{0}\left(\epsilon_{{\rm tot}}+e\phi({\bf r})\right)$
(band bending approximation) and the bottom of the conduction band
now corresponds to an energy $E_{C}-e\phi({\bf r})$. The electrochemical
potential $\bar{\mu}=\mu_{0}({\bf r})-e\phi({\bf r})$ is uniform
through the system, and the chemical potential thus varies as $\mu_{0}({\bf r})=\bar{\mu}+e\phi({\bf r})$.
The local density of particles depends on the (uniform) electrochemical
potential and the local density of state:
\begin{align*}
n({\bf r}) & =\int d\epsilon_{{\rm tot}}\,\frac{D(\epsilon_{{\rm tot}},{\bf r})}{\exp\left(\frac{\epsilon_{{\rm tot}}-\bar{\mu}}{k_{B}T}\right)-1}
\end{align*}
Note that we can also express the density of carriers as a function
of the local chemical potential $\mu_{0}({\bf r})$ and the unshifted
(uniform) density of state:
\begin{equation}
n({\bf r})=\int d\epsilon\,\frac{D_{0}(\epsilon)}{\exp\left(\frac{\epsilon-\mu({\bf r})}{k_{B}T}\right)-1}
\end{equation}

A: Disclaimer : I'm not sure that the following is the exact/complete answer to the question but maybe these elements could help.

Equilibrium of a system under an external field

Let say that you have an open macroscopic system $\Sigma$ (composed of identical particles) which is under the influence of an external time-independant but space-dependant potentiel field $\phi(\textbf{x})$. One can slice the entire system in an ensemble of subsytems $\{\delta\Sigma\}$. At this point, $\phi$ is assumed to be smooth enough function of $\textbf{x}$ so that it can be considered as constant at the scale of these $\delta\Sigma$'s.
In the absence of external force, the particles contained in each $\delta\Sigma$ are described by the hamiltonian $\hat{\mathcal{H}}_0$ which eigenstates $\{\vert\ell\rangle\}_\ell$ are supposed to be known, giving then the corresponding eigen energies $E_\ell$. These $\vert\ell\rangle$ states should be eigenstates of the particle number operator $\hat{N}$, associated to the eigenvalue $N_\ell$. In such case, the grand canonical partition function reads :
$$
\Xi_0(\beta,\mu)=\text{Tr}\,e^{-\beta(\hat{\mathcal{H}}_0-\mu\hat{N})}
$$
where $\beta$ is the temperature, $\mu$ is the chemical potential.
The equilibrium condition states here that the temperature $\beta$ and $\mu$ have to be constant in the all system (the same in each $\delta\Sigma$).
Now in the presence of an external potential $\phi(\textbf{x})$, the new hamiltonian is $\hat{\mathcal{H}}=\hat{\mathcal{H}}_0+\hat{N}\phi(\textbf{x})$. For each sub-system $\delta\Sigma$, each $\vert\ell\rangle$ state has now an associated eigen energy
$$
E'_\ell=E_\ell+N_\ell\phi(\textbf{x})
$$
In this case the partition function is then :
$$
\Xi(\beta,\mu,\textbf{x})=\text{Tr}\,e^{-\beta(\hat{\mathcal{H}}_0-\bar{\mu}(\textbf{x})\hat{N})}\equiv\Xi_0(\beta,\bar{\mu}(\textbf{x}))
$$
where $\bar{\mu}(\textbf{x})=\mu-\phi(\textbf{x})$ is the local chemical potential. In that case, one can show that the equilibrium condition still constant chemical potential $\mu$ and temperature $\beta$. In particular, this implies that :
$$
\forall\textbf{x},\,\bar{\mu}(\textbf{x})+\phi(\textbf{x})=\mu=c^{st}
$$

Application for an electronic gas under an external electric field $\textbf{E}_{ext}$

In this case, the external electrostatic potential $V_{ext}$ (defined as $\textbf{E}_{ext}=-\nabla V_{ext}$) is associated to a local charge distribution $\rho_{ext}$ through a Poisson equation :
$$
\Delta V_{ext}+\frac{\rho_{ext}}{\epsilon_0}=0
$$
In reaction to this external potential, electron will generate an internal charge density $\rho_{in}$ (this mechanism is at the origin of the screening effects). To $\rho_{in}$ can be associated an electrostatic potential $V_{in}$ via a Poisson equation as well. 
Finally, the total electrostatic potential inside the metal is :
$$
V(\textbf{x})=V_{ext}(\textbf{x})+V_{in}(\textbf{x})
$$
The equilibrium condition is then :
$$
\forall\textbf{x},\,\bar{\mu}(\textbf{x})-eV(\textbf{x})=\mu=c^{st}
$$

The origin of the confusion?

I think that this last equation maybe the origin of the confusion. A priori, in the way it is defined in the partition function $\Xi$, the local chemical potential $\bar{\mu}$ does entail an electrostatic contribution. Conversly, by definition $\mu$ does not entail any electrostatic contribution.
However, in the special case where one considere an equilibrium configuration, the only intensive experimental accessible value of the chemical potential is $\mu$ (not $\bar{\mu}$). And in that case, spatial variations of $\bar{\mu}$ exactly compensate variations of $V$ in order to have a constant $\mu$. In that very special case I think that people would call $\mu$ the "electrochemical potential" because of the $-eV(\textbf{x})$ contribution in the last equation. But it should be kept in mind that it is just an artifact of considering an equilibrium state of the system.
So now, in the background of semiconductors physics, the following has quite convinced myself that the Fermi level $E_F$ is $\mu$ and not $\bar{\mu}$.
