With suitable boundary conditions, the decomposition is unique. Without them, it's not.
Suppose that $(\phi,{\bf G})$ and $(\phi',{\bf G}')$ are two different decompositions for the same function. Then
$$
\nabla(\phi-\phi')+\nabla\times({\bf G}-{\bf G}')=0.
$$
Take the divergence of both sides to find that
$$
\nabla^2(\phi-\phi')=0.
$$
So for any two distinct decompositions, the scalar field $\phi$ must differ by a harmonic function $f$ (that is, one with $\nabla^2f=0$). Moreover, any harmonic function will work -- that is, there will be a way to choose a ${\bf G}'$ to go along with this $\phi'$. To see this, note that we have to choose ${\bf G}'$ to satisfy
$$
\nabla\times({\bf G}'-{\bf G})=\nabla f.
$$
The right side of this expression is divergence-free (because it's $\nabla^2f$), and any divergence-free vector field can be expressed as the curl of some other divergence-free vector field, so ${\bf G'}-{\bf G}$ exists.
(A couple of notes: This latter fact is the one that lets us define the vector potential for a given magnetic field, specifically in Coulomb gauge. To be honest, I don't remember the proof that there exists a function ${\bf G}$ whose curl is ${\bf B}$ for any divergence-free ${\bf B}$. I do remember how you show that, having gotten such a ${\bf G}$, you can make it divergence-free: Just subtract off $\nabla q$ where $\nabla^2q=\nabla\cdot{\bf G}$. The new ${\bf G}$ will have the same curl as the old one and will be divergence-free.
One other thing: complications arise if the domain we're considering isn't simply connected. Let's say it is.)
So the answer is that, to make the decomposition unique, you have to impose strong enough boundary conditions to make it so that no harmonic functions exist. For a compact domain without boundary (such as the surface of a sphere), you don't need any boundary conditions: there are no non-constant harmonic functions on such domains. (Slick proof of this: you can prove that harmonic functions never have local maxima or minima, but a nonconstant function on such a domain must have them -- in particular, it must have a global maximum and a global minimum somewhere.)
For a compact region with boundary, you need to specify either $\phi$ or the normal component of $\nabla\phi$ on the boundary. For good old infinite space, you need to specify that $\phi$ approach zero (or some other given function) as you tend to infinite distance.
It's easy to check that without such boundary conditions, you get into trouble. For instance, take the functions
$$
\phi=x, {\bf G}=z\hat{\bf j}.
$$
They give rise to
$$
{\bf H}=\nabla\phi+\nabla\times{\bf G}=\hat{\bf i}-\hat{\bf i}=0.
$$
So this pair can be added to any Helmholtz decomposition without changing the original vector field.
Let me try this more clearly than the other answers, which aren't wrong. You ask:
So, can someone please elaborate what this EM field is with respect to $\vec E$ and $\vec B$ in the context of Helmholtz decomposition?
There is no "EM field in the context of Helmholtz decomposition". Helmholtz just says that every vector field $\vec V$ is decomposable as curl and gradient of two other fields, i.e.
$$\vec V = \vec \nabla \phi + \vec \nabla \times \vec A $$
You can do this for the electric or the magnetic field, of course, but this isn't particularly enlightening as to the nature of "the EM field". A field should behave nicely under transformations, and special relativity with its action on the electric and magnetic fields shows us that we should not add them together, but seek a quantity that transforms nicely under Lorentz transformations instead:
"The electromagnetic field" is equivalently the gauge four-potential $A$ (consisting of the scalar electrostatic potential in the temporal and the magnetic vector potential in the spatial entries) or its derivative, the field strength tensor $F = \mathrm{d}A$. Electric and magnetic fields become part of the tensor as
\begin{align}
F^{0i} & = E^i \\
F^{ij} & = \sum_k\epsilon^{ijk}B^k
\end{align}
This is "the EM field", but it has nothing to do with Helmholtz decomposition, since electromagnetism is properly looked at in the four-dimensional setting of special relativity, for which only the general Hodge decomposition may be applied, of which Helmholtz is a special case, but even this has nothing to do with it.
This EM field acts on the four-velocity, reproducing the Lorentz force by
$$ \frac{\mathrm{d}p}{\mathrm{d}t} = q F(u)$$
where $u$ is the four-velocity, and $(F(u))_\mu = F_{\mu\nu}u^\nu$.
Best Answer
This is answers follows from the comments to the initial text. First of all I think there is one major problem with the question:
It is too general. Helmholtz Decomposition (HD) is a mathematical technique. There is no need to use it, in principle, but often it can make life easier. The best way to learn about the cases where HD is useful is to do specific examples, rather than ask general questions. For example, there are books by Griffiths and Jackson (Classical Electrodynamics). There are also more rigorous treatments by Morse & Feshbach for example. Going towards even more maths and less physics, the Helmholtz theorem is treated by Arfken & Weber (Mathematical Methods...).
Having said this, I can try to give my best personal view why HD is useful. Lets us take the statement of the theorem by Arfken & Weber (Ch3. Helmholtz's Theorem).
* A vector field is uniquely specified by giving its divergence and its curl within a simply connected region and its normal component on the boundary. *
So let us look at the electrostatic Maxwell's Equations.
$$\boldsymbol{\nabla}\cdot \mathbf{E}=\frac{\rho}{\epsilon}$$
$$\boldsymbol{\nabla}\times\mathbf{E}=\mathbf{0}$$
Helmholtz theorem tells us that these equations are sufficient to fully and uniquely constrain the solution for electric field ($\mathbf{E}$) given boundary conditions. So any solution you can find for this problem is the only correct one.
Better still the Wikipedia link I gave (https://en.wikipedia.org/wiki/Helmholtz_decomposition) goes further and states that electric field can be written as:
$$\mathbf{E}=-\boldsymbol{\nabla}\phi + \boldsymbol{\nabla}\times\mathbf{K}$$
and then:
$$\phi\left(\mathbf{r}\right) = \frac{1}{4\pi}\int_V d^3 r' \frac{\boldsymbol{\nabla'}\cdot\mathbf{E}\left(\mathbf{r}'\right)}{\left|\mathbf{r}-\mathbf{r}'\right|} - \frac{1}{4\pi}\oint_{\partial V} d^2 r'\frac{\mathbf{\hat{n}'}\cdot\mathbf{E}\left(\mathbf{r}'\right)}{\left|\mathbf{r}-\mathbf{r}'\right|}$$
Where $V$ is the simply connected region. Provided $\mathbf{E}$ vanishes at the surface of $V$ (i.e. at $\partial V$) faster than $1/r$ (here is your boundary condition), we get:
$$\phi\left(\mathbf{r}\right) = \frac{1}{4\pi}\int_V d^3 r' \frac{\rho\left(\mathbf{r}'\right)/\epsilon}{\left|\mathbf{r}-\mathbf{r}'\right|} $$
Similar treatment exists for $\mathbf{K}$, see link. So using Helmholtz decomposition you can find the electrostatic field.