Both methods are "appropriate".
IDW basically says "half way between a point valued 6 and a point valued 8 the value is going to be 7", but extended to the distance to multiple points. It can't say anything about how confident it is that the value is 7 though. It gives a nicely smooth surface through the data.
Kriging basically says "your data come from a correlated surface and I'll work out the correlation from pairs of values, and then work out the mean at each interpolated location and use the correlation to also work out how unsure I am about that mean. So I think its around 7, and between 6.5 and 7.5 with 95% confidence".
Your spaced points will make it hard for Kriging to work out the short-range correlation, so its estimates will probably head towards the global mean value close to your sample points, with a large uncertainty. Kriging is honest, its saying "I don't know what the value is 100m from your data point, but all your data points average to 18.06 with a variance of 0.04, so... maybe that?" Nearer to sample points it will have a smaller variance, and at the sample point it will be zero because it knows what the value was there.
Another option is to fit a polynomial in x and y to some number of degrees (Z = A + Bx + Cy + Dxy + Ex^2 + Fy^2 + ...) and that will give you a model object you can evaluate over the space to get a smooth curve, but it will go pretty crazy outside the space (quadratically) so crop it tightly to the region.
But it all comes down to why you want to do this, and how concerned you are about accuracy and precision.
I've just worked through very simple applications of kriging, IDW and polynomial surface fitting and got these:
Kriging:
As expected, there's not much in the variogram to assess the spatial correlation, so kriging says "I don't really know what it is anywhere not near any of the points, so my estimate is 'its the average' and my standard error is 'large'".
IDW:
Very glamorous and seductive, in that it looks like a plausible underlying surface, but in reality all sorts of noise could really be in there and it gives no estimate of the uncertainty.
Polynomial smoothing surface:
With only nine points you can't put a lot of wigglyness in a polynomial surface, you don't have that many degrees of freedom to play with. This fitted surface doesn't go through the points, being a smoothing rather than an interpolation.
R code follows. I know you tagged ESRI tools but this is more about the method comparison than tool application so consider these an appendix. ESRI tools will probably do similar outputs.
I did a sort of fake georeferencing of the data and digitised it, so this is approximate. Saved as a geopackage and read in with sf but then converted to sp because the software needed that. Required some faffing with the CRS to make it work. Couldn't fit a y^2 without an NA.
library(sf)
library(sp)
library(ggplot2)
library(automap)
library(gstat)
q = st_read("./q.gpkg")
qsp = as(q,"Spatial")
proj4string(qsp) = "+init=epsg:27700"
comment(proj4string(qsp)) <- NULL
## Kriging
k = autoKrige(ZZ~1, input_data=qsp)
plot(k)
## IDW
id = idw(ZZ~1, qsp, newdata=k$krige_output)
plot(id)
## Model
qpp = cbind.data.frame(coordinates(qsp), ZZ=qsp$ZZ)
names(qpp)[1:2] = c("x","y")
m0 = lm(ZZ ~ x+y+I(x^2)+I(x*y), data=qpp)
newdata = data.frame(coordinates(k$krige_output))
names(newdata) = c("x","y")
newdata$fitted = predict(m0, newdata)
ggplot(newdata) + geom_point(aes(x=x,y=y,col=fitted))
and here's the data in CSV. Left as an exercise to read this in and convert to geopackage or for you to work with it directly and convert to sp spatial data frame:
"x","y","ZZ"
348063.028597848,458037.963664607,17.910456
348290.770878634,458027.026261644,17.922374
348421.122149715,458007.82288224,18.137663
348051.993373532,457841.866673889,17.961751
348237.073696967,457785.449070263,18.009223
348410.623309634,457831.170723402,18.179926
348037.300667433,457661.972826257,17.828845
348204.641407268,457631.404998955,17.756329
348383.262432327,457628.681298612,17.897605
