As long as you're producing the surface from a matrix of values, you could use plotly to do something like this:
# 2D kernel density estimation
kd <- with(MASS::geyser, MASS::kde2d(duration, waiting, n = 50))
x <- kd$x
y <- kd$y
# this was rows correspond to y and columns corresponds to x
# (this is the way plotly prefers)
z <- t(kd$z)
plot_ly(x = x, y = y, z = z, type = "surface", size = I(5)) %>%
# condition on the 10th x value
add_paths(x = x[10], y = y, z = z[, 10], color = I("red"), linetype = I("dash")) %>%
# project it to the 'background'
add_paths(x = min(x), y = y, z = z[, 10], color = I("black")) %>%
# condition on the 15th y value
add_paths(x = x, y = y[15], z = z[15, ], color = I("red"), linetype = I("dash")) %>%
# project it to the 'background'
add_paths(x = x, y = min(y), z = z[15, ], color = I("black"))