  In this paper, we continue to consider the generalized Liouville system: $$
\Delta_g u_i+\sum_{j=1}^n a_{ij}\rho_j\left(\frac{h_j e^{u_j}}{\int h_j
e^{u_j}}- {1} \right)=0\quad\text{in \,}M,\quad i\in I=\{1,\cdots,n\}, $$ where
$(M,g)$ is a Riemann surface $M$ with volume $1$, $h_1,..,h_n$ are positive
smooth functions and $\rho_j\in \mathbb R^+$($j\in I$). In previous works
Lin-Zhang identified a family of hyper-surfaces $\Gamma_N$ and proved a priori
estimates for $\rho=(\rho_1,..,\rho_n)$ in areas separated by $\Gamma_N$. Later
Lin-Zhang also calculated the leading term of $\rho^k-\rho$ where $\rho\in
\Gamma_1$ is the limit of $\rho^k$ on $\Gamma_1$ and $\rho^k$ is the parameter
of a bubbling sequence. This leading term is particularly important for
applications but it is very hard to be identified if $\rho^k$ tends to a higher
order hypersurface $\Gamma_N$ ($N>1$). Over the years numerous attempts have
failed but in this article we overcome all the stumbling blocks and completely
solve the problem under the most general context: We not only capture the
leading terms of $\rho^k-\rho\in \Gamma_N$, but also reveal new robustness
relations of coefficient functions at different blowup points.
