Blow-up solutions for mean field equations with Neumann boundary conditions on Riemann surfaces

On a compact Riemann surface $(\mathrm{\Sigma },g)$ with a smooth boundary $\mathrm{\partial }\mathrm{\Sigma }$, we consider the following mean field equations with Neumann boundary conditions: $$-{\mathrm{\Delta }}_{g}u=\lambda (\frac{V{e}^u}{{\int }_{\mathrm{\Sigma }}V{e}^{u}d{v}_g}-\frac{1}{\mathrm{\mid }\mathrm{\Sigma }{\mathrm{\mid }}_g})\text{ in }\mathrm{\Sigma }\text{ with }{\mathrm{\partial }}_{{\nu }_g}u=0\text{ on }\mathrm{\partial }\mathrm{\Sigma },$$ We find conditions on the potential function $V:\mathrm{\Sigma }\to {\mathbb{R}}^{+}$ such that solutions exist for the parameter $\lambda$ when it is in a small right (or left) neighborhood of a critical value $4\pi (m+k)$ for $k\le m\in {\mathbb{N}}_{+}$ and blow up as $\lambda$ approaches the critical parameter. The blow-up occurs exactly at $k$ points in the interior of $\mathrm{\Sigma }$ and $(m-k)$ points on the boundary $\mathrm{\partial }\mathrm{\Sigma }$.