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.
Your observation is correct. For steady irrotational flow of an incompressible, inviscid fluid with constant density, the Euler equation of momentum reduces to
$$\nabla \left(\frac{p}{\rho} + \frac{1}{2} |\mathbf{u}|^2 + \chi \right) = 0$$
Integrating the components with respect to the spatial variables, we get the general solution
$$\frac{p}{\rho} + \frac{1}{2} |\mathbf{u}|^2 + \chi = c(t),$$
where the arbitrary function $t \mapsto c(t)$ changes nothing about the flow field and can be absorbed into the pressure. Recall that the pressure field corresponding to a particular velocity field is never uniquely determined since an arbitrary reference pressure can be added. The pressure and body-force potential only affect the velocity through their spatial gradients and adding an arbitrary function of time has no impact.
More generally for irrotational flow where $\nabla \times \mathbf{u} =0$, the velocity field is the gradient of a potential, $\mathbf{u} = \nabla\phi$, and the momentum equation reduces to
$$\frac{\partial\nabla \phi}{\partial t}= -\nabla \left(\frac{p}{\rho} + \frac{1}{2} |\nabla \phi|^2 + \chi \right),$$
which upon integration yields
$$\frac{\partial \phi}{\partial t}+\frac{p}{\rho} + \frac{1}{2} |\nabla \phi|^2 + \chi = c(t)$$
Again, in steady flow where $\frac{\partial \phi}{\partial t} = 0$, the appearance of the time-dependent $c(t)$ is not excluded.
Best Answer
The key concept needed here is that the Hemholtz decomposition is not necessarily unique.
Non-uniqueness can occur because there exist nontrivial vector fields which are both irrotational and divergence-free. For example, the constant 2D velocity field $\vec v = (1,0)$ can be expressed as either $\vec v=-\nabla \phi$ with $\phi(x,y)=-x$, or as $\vec v=\nabla \times \vec \Psi$ with $\vec \Psi (x,y) = (0,0,y)$.
Because of this non-uniqueness, showing that it’s possible to express a particular velocity field with an incompressible flow as a Hemholtz decomposition with a nonzero $\phi$ doesn’t mean that it isn’t also possible to express the same velocity field as a different Hemholtz decomposition in which $\phi = 0$.
An arbitrary 2-D velocity field with $\nabla \cdot \vec v = 0$ can be written purely in terms of a vector potential in which the stream function is
$$\psi(x,y)=\int^{y}_{0}v_{x}(0,y')dy'-\int^{x}_{0}v_{y}(x',y)dx' .$$
I'll leave verification of that equation as an exercise to the reader.