In questions that ask about Planck's constant entering into statistical mechanics, a common and accepted answer is that Planck's constant is an arbitrary normalization that falls out when calculating experimentally measurable quantities.
Specifically, it's said in the above questions that $h$ enters to de-dimensionalize the product $dp dq$ or to normalize $dp dq$ to count states, and the choice for such a constant with units of action is nearly arbitrary.
To the contrary, consider how one calculates quantities in the grand canonical ensemble with a variable number of particles - mentioned in this question.
In the grand canonical ensemble, the key potential is the grand potential, which is $$\Phi = U - TS-\mu N .$$ The link to statistical mechanics comes from $$\Phi = -kT \ln(\mathcal{Z}), $$
where $$\mathcal{Z} = \sum_N e^{\beta \mu N} Z(V,N,T).$$
For simplicity, consider a classical partition function $Z$ for N non-interacting particles. We have then that
$$Z = \frac{1}{h^{3N} N!} \Big(\int d^3p d^3q \, e^{-\beta E({p}, {q})}\Big)^N$$
This yields in turn that
$$\mathcal{Z} = e^{\frac{1}{h^3} \Big(\int d^3p d^3q \, e^{-\beta E({p}, {q})}\Big)e^{\beta \mu}}$$
and in turn that $\Phi$ is proportional to $\frac{1}{h^3}$:
$$\Phi = -kT \frac{1}{h^3} \Big(\int d^3p d^3q \, e^{- \frac{E({p}, {q})}{kT}}\Big)e^{\frac{\mu}{kT}}. $$
For example, for non-interacting massless particles with $E = cp$ in a box with two internal degrees of freedom with $\mu = 0$ (i.e. a naive picture of light, without any $h$ put into the energy or whatnot), we have
$$\Phi = -\frac{16 \pi k^4}{h^3 c^3} T^4 V.$$ This leads to a pressure
$$P = \frac{16 \pi k^4}{h^3 c^3} T^4$$ from the thermodynamic relation above.
As you can see, $h$ enters into an experimentally measurable quantity! Thus one could experimentally measure Planck's constant via the pressure of such a gas. I'll leave it to you to check that $h$ also enters into the pressure of a non-relativistic gas at a fixed $\mu$.
I will take from the above arguments that $h$ is not simply an arbitrary constant for de-dimensionalizing the product $dp dq$. Given this point of view that the de-dimensionalization of $dp dq$ has experimental ramifications, can we show that for every system in thermodynamic equilibrium that the de-dimensionalizing constant $h$ must be a universal constant? That is, since the de-dimensionalization is not arbitrary, can we show the de-dimensionalizing constant must be the same for all equilibrium systems in classical statistical mechanics?