g[x_,y_]=Exp[-(x^2)-(y^2)]

{D[g[x,y],x],D[g[x,y],y],-1}

n=%/.{x->0,y->0}

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

{d1,d2}={{1,0,D[g[x,y],x]},{0,1,D[g[x,y],y]}}/.{x->0,y->0}

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

Show[b,c]

n=Cross[d1,d2];

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

Show[b,c,d]