a=ContourPlot3D[z == x^2 + y^2, {x,-3,3}, {y,-3,3},{z,0,9},ContourStyle->Opacity[0.2],Boxed->False]

b=Show[ContourPlot3D[z == x^2 + y^2, {x,-3,3}, {y,-3,3},{z,0,9},ContourStyle->Opacity[0.2],Boxed->False, MeshFunctions->{#1&,#2&},Mesh->{{{2}},{{1}}},MeshStyle->{{Thick,Blue},{Thick,Red}}],Graphics3D[{PointSize[Large],Black,Thick,Point[{2,1,5}]}],ContourPlot3D[{x==2,y==1},{x,-3,3}, {y,-3,3},{z,0,9},Mesh->False,ContourStyle->{Opacity[0.1],Opacity[0.1]}]]

g[x_,y_]=x^2+y^2;

Grad[g[x,y],{x,y}];

{d11,d21}={{1,0,%[[1]]},{0,1,%[[2]]}};

%/.{x->2,y->1}

{d1,d2}={(1/Norm[%[[1]]])%[[1]],(1/Norm[%[[2]]])%[[2]]}

c=ParametricPlot3D[{{2,1,5}+t d1,{2,1,5}+t d2},{t,-2,2},PlotStyle->{{Thick,Blue},{Thick,Red}}];

Show[b,c]

n=Cross[d1,d2];

d=ContourPlot3D[n.{x,y,z}==n.{2,1,5}, {x,1,3}, {y,0,2},{z,4,6},Mesh->False,ContourStyle->{LightGray,Opacity[0.4]},BoundaryStyle->{Black,Thick}];

Show[b,c,d]