Your understanding is correct; you need to have a well-defined reference level for the Fermi level if you want to make comparisons between different systems (for example, comparing work functions of different surfaces for the same metal or different metals). In addition to subtle variations due to pseudopotential selection, you'll also have different Fermi level / work function predictions depending on how many layers of material you include, how many layers you are freezing, whether you are using a symmetric or asymmetric slab, etc. For example, for metals, it's well-known that the work function converges somewhat slowly in an oscillatory fashion with respect to layer thickness.
For condensed phases with 3D periodicity, you don't have a well-defined reference level. You also cannot have a bulk system that has a finite charge; it's electrostatic potential energy is undefined (in reciprocal space, the electrostatic potential diverges at the gamma point) since you'd have a periodic collection of charged cells. Most planewave DFT codes will insert a uniform background charge (equivalent to setting the electrostatic potential at the gamma point to zero) to cancel the net charge to avoid this and still enable the calculation to converge. For these reasons, you couldn't build an electrochemical model for a bulk (3D periodic) structure or easily use a GC-SCF method.
But essentially, the work function and electrode potential are inherently surface-dependent properties because they require the presence of an interface to even define them in the first place. You'd also need to include some sort of compensating charge distribution (like an electrical double layer) to keep the system at a net neutral charge.