It is correct that equivalence of ensemble fails in general at a phase transition (but see the discussion at the bottom of this answer). However, phase separation can perfectly well be described using statistical mechanics, even mathematically rigorously. Let me briefly explain this, in the simpler setting of an Ising lattice gas. (Here is a review on this topic.)
Consider an Ising model in the box $\Lambda_N=\{1,\dots,N\}^d$, with $+$ boundary condition. Let $m_N = \frac{1}{N^d}\sum_{i\in\Lambda_N} \sigma_i$ be the magnetization in the box. Consider the Gibbs measure conditioned on the event $m_N = m \in (-m^*(\beta),m^*(\beta))$, where $m^*(\beta)$ is the spontaneous magnetization at inverse temperature $\beta$. The constraint can be seen as considering the Ising lattice gas (say, mapping $-$ spins to particles and $+$ spins to holes) in the canonical ensemble, that is, with a fixed number of particles in the box (equal to $\frac{1-m}{2}N^d$).
Then, one can prove that the following event occurs with a probability converging to $1$ as $N\to\infty$:
- There is a single giant droplet of $-$ phase immersed inside the $+$ phase.
- Once scaled by $1/N$, the shape of this droplet converges (as $N\to\infty$) to the deterministic shape that minimizes the surface tension and of the volume required to satisfy the constraint. (More precisely,the latter set is the solution to the following variational problem: it is the set $V\subset [0,1]^d$ with volume $(m^*(\beta)-m)/2m^*(\beta)$ that minimizes $\int_{\partial V} \tau(\vec n_s) \mathrm{d}\mathcal{H}^{d-1}_s$, where $\tau(\vec n)$ is the surface tension in direction $\vec n$, $\vec n_s$ is the normal to $V$ at $s$ and $\mathcal{H}^{d-1}$ is the $(d-1)$-dimensional Hausdorff measure.)
A picture being worth a thousand words, here is a typical configuration of the 2d model (particles, i.e. $-$ spins, are in dark blue):
Concerning your questions
However, at zero Celsius of water-ice mixture, is boltzmann ensemble
not equivalent to microcanonical ensemble?
Is the Boltzmann ensemble not sharp in energy per particle but is an
integration of microcanonical ensemble over a big range of energy?
Note that the energy densities of both phases are equal in the above model, so that the only effect of the phase separation is to introduce a surface order correction to the energy. The latter, of course does not impact the energy per particle (in a large system). I'd guess, although I haven't really thought hard about it, that equivalence of the above canonical ensemble and the corresponding microcanonical ensemble should therefore hold in this case. This has never been proved in this setting, as far as I know.
The lack of equivalence of ensembles mentioned at the beginning of this answer would be seen, in this case, between the canonical and the grand canonical ensembles: for all values of the imposed magnetization $m\in (-m^*(\beta),m^*(\beta))$ in the canonical ensemble, the conjugate magnetic field $h$ is equal to $0$ and typical configurations in the grand canonical ensemble (at $h=0$) would not display the type of phase separation described above (what would happen depends on the choice of boundary condition; for the $+$ boundary condition chosen above, typical configurations would correspond to a homogeneous $+$ phase).
