Bardeen-Zumino polynomials

Question

Suppose the theory with chiral fermions $\psi$ and gauge fields $B$. In general, it contains the chiral anomaly, due to which the VEV $\langle J_{\mu}(x)\rangle_{B}$ has anomalous divergence. When we talk about how to define the anomaly $\text{A}[B]$, there are two different ways. Suppose the quantum effective action $\Gamma[B]$,

$$
e^{i\Gamma[B]} = \int D\bar{\psi}D\psi \text{exp}\left[ i\int\bar{\psi}i\gamma_{\mu}D^{\mu}\psi - \frac{1}{4}F_{\mu\nu}^{a}F^{\mu\nu,a}\right]
​$$

Then the consistent anomaly is defined such that

$$
\tag 1 \delta_{\epsilon}\Gamma [B] = \text{A}[B], \quad \langle J_{\mu,b}(x)\rangle^{\text{consistent}}_{B} \equiv \frac{\delta \Gamma[B]}{\delta B^{\mu}_{b}},
​$$

so that

$$
\partial_{\mu}\langle J^{\mu}_{b}\rangle_{B}^{\text{consistent}} = \text{A}^{\text{consistent}}_{b}
​$$

It states that the anomaly is defined through the gauge variation of the quantum effective action; therefore the consistent anomaly automatically satisfies the Wess-Zumino consistency condition.

On the other hand, it can be shown that the resulting anomaly $\text{A}^{\text{consistent}}(x)$ isn't gauge covariant, which states also that a covariant divergence of the current isn't gauge covariant. If this current is gauge one, it is preferable to choose it in a way such that its divergence is gauge covariant.  This can be done by modifying $(1)$ by changing the definition of current:

$$
\langle J_{\mu,b}\rangle^{\text{covariant}}_{B} \equiv \frac{\delta \Gamma[B]}{\delta B^{\mu}_{b}} + J_{BZ},
​$$

where $J_{BZ}$ is the so-called Bardeen-Zumino polynomial - the local gauge-variant expression of gauge fields, which leads to the modification of the anomalous conservation law, making it covariant:
$$
\partial_{\mu}\langle J^{\mu}_{b}(x)\rangle^{\text{covariant}}_{B} = \text{A}^{\text{covariant}}_{b}(x)
​$$

My question is: what is the reason for ambiguity of the definition of the current so that we can add the Bardeen-Zumino polynomial? It's not the regularization procedure arbitrariness, if I understand correctly, because this arbitrariness only provide us to choose the coefficient in the front of the anomaly in cases when there are different currents in triangular vertex; also the counterterm in the effective action, which regulates the coefficient in the front of the anomaly, has to be the local function of the gauge fields, while the BZ polynomial can't be obtained as the gauge variation of effective action.