Daniel,
First,  to get your bar legend, you need to make a change in syntax:
Plot3D[Re[fun^2 + 1], {x, -2, 2}, {y, -2, 2}, 
 ColorFunction -> Function[{x, y, z}, Hue[Im[(x + y I)^2 + 1]]], 
 PlotLegends -> BarLegend[{Hue, {-4, 4}}]]
As to the Color function, the problem is the default behavior is ColorFunctionScaling->True which scales your x,y, and z to the range 0 to 1.  This means that your calculation of Im[fun^2+1] is no longer correct. You need to use ColorFunctionScaling->False, and then you need to scale the values of the argument to Hue to make sense (between 0 and 1). This means you need an offset and scale to Hue to get your function values between 0 and 1.  You can automate this easily but the basic plot should look something like this:
Plot3D[Re[fun^2 + 1], {x, -2, 2}, {y, -2, 2}, 
 ColorFunction -> 
  Function[{x, y, z}, Hue[(Im[(x + y*I)^2 + 1] + 8)/16]], 
 ColorFunctionScaling -> False, PlotLegends -> BarLegend[Hue]]
The legend has the wrong values but you can fix that by creating a function to give to BarLegend that scales the values so the bar is correct.
You get the following:

Regards,
Neil