Shape of a barchan dune
reference : G. Sauermann, P. Rognon, A. Poliakov, H.J. Herrmann : "The shape of the barchan duns of Southern Morocco", Geomorphology 36, 47 - 62 (2000)  sauermann_barcan_shape.pdf

all dimensions in meters

implement instructions in red and cyan background

Define the principal functions of dune height using constants found in the paper above :

l[h_]:=al*h+bl (* dune length upwind from origin al = 5.4  bl = 10.7*);
w[h_]:=aw*h+bw (* dune width aw = 5.8  bw = 1.9 *);
b0[h_]:=ab0*h+bb0 (* brink curve coefficient ab0 = -2.4  bb0 = 7.9 *);
b2[h_]:=bb2*l[h]/w[h]^2 (* brink curve coefficient bb2 = -1.46 *);
zz[x_,y_,h_]:=h*(1-(x/l[h])^2-(y/w[h])^2) (* upwind dune height *);
bx[y_,h_]:=b0[h]+b2[h]*y^2 (* x-projection of the brink curve on the ground *);
bz[y_,h_]:=zz[bx[y,h],y,h] (* height of the brink curve *);

Simplify[(∂_y) bx[y,h]]

derbx[y_,h_]:=((2 bb2 (bl+al h) y)/(((bw+aw h))^2)) (* derivative the brink curve *);

bx0[y_,z_,h_]:=bx[y,h]-(bz[y,h]-z)/Tan[alpha]/Sqrt[1+derbx[y,h]^2] (* x-intersection of the slip face with height z starting from a location of y at the brink; alpha = Pi/4 typically *);
by0[y_,z_,h_]:=y+(bz[y,h]-z)/Tan[alpha]*derbx[y,h]/Sqrt[1+derbx[y,h]^2] (* y-intersection of the slip face with height z *);

Calculate the intersection of the brink and the windward face :

Solve[(bx[y,h]/l[h])^2+(y/w[h])^2==1,y]

Simplify[√(-((bb0 bb2 bl (bw^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((bl^2) (bw^2))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2))))-((2 aw bb0 bb2 bl bw h)/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((aw (bl^2) bw h)/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((al bb0 bb2 (bw^2) h)/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((al bl (bw^2) h)/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((ab0 bb2 bl (bw^2) h)/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((aw^2) bb0 bb2 bl (h^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((aw^2) (bl^2) (h^2))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2))))-((2 al aw bb0 bb2 bw (h^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((2 al aw bl bw (h^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((2 ab0 aw bb2 bl bw (h^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((al^2) (bw^2) (h^2))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2))))-((ab0 al bb2 (bw^2) (h^2))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((al (aw^2) bb0 bb2 (h^3))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((al (aw^2) bl (h^3))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((ab0 (aw^2) bb2 bl (h^3))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((al^2) aw bw (h^3))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-((2 ab0 al aw bb2 bw (h^3))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))-(((al^2) (aw^2) (h^4))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2))))-((ab0 al (aw^2) bb2 (h^4))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))+(((bw^2) (((bl+al h))^(3/2)) (Sqrt[4 bb0 bb2+bl+4 (bb2^2) bl+al h+4 ab0 bb2 h+4 al (bb2^2) h]))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2))))+((aw bw h (((bl+al h))^(3/2)) (Sqrt[4 bb0 bb2+bl+4 (bb2^2) bl+al h+4 ab0 bb2 h+4 al (bb2^2) h]))/((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))+(((aw^2) (h^2) (((bl+al h))^(3/2)) (Sqrt[4 bb0 bb2+bl+4 (bb2^2) bl+al h+4 ab0 bb2 h+4 al (bb2^2) h]))/(2 ((bb2^2) (bl^2)+2 al (bb2^2) bl h+(al^2) (bb2^2) (h^2)))))]

ym4[h_]:=(Sqrt[-(((((bw+aw h))^2) (2 bb0 bb2+bl+al h+2 ab0 bb2 h-(Sqrt[bl+al h]) (Sqrt[4 bb0 bb2+bl+4 (bb2^2) bl+(al+4 ab0 bb2+4 al (bb2^2)) h])))/(2*(bb2^2) (bl+al h)))]) (* positive y-coor of the outermost point of the slip face on the ground = intersection of the brink and the windward face *);
xm[h_]:=bx[ym4[h],h] (* x-coor of that point *);

Simplify[Series[xm[h],{h,0,2}]]

Solve for height on the slip face of a given point at (x, y):

Solve[yy-y+derbx[y,h]*(xx-bx[y,h])==0,y] (* find the y-coord at the brink above the point of interest at (xx,y) *)

ybrink[xx_,yy_,h_]:=Re[-(2 (6^(1/3)) bb0 (bb2^3) (bl^3) (bw^2)+(6^(1/3)) (bb2^2) (bl^2) (bw^4)+4 (6^(1/3)) aw bb0 (bb2^3) (bl^3) bw h+6 (6^(1/3)) al bb0 (bb2^3) (bl^2) (bw^2) h+2 (6^(1/3)) ab0 (bb2^3) (bl^3) (bw^2) h+4 (6^(1/3)) aw (bb2^2) (bl^2) (bw^3) h+2 (6^(1/3)) al (bb2^2) bl (bw^4) h+2 (6^(1/3)) (aw^2) bb0 (bb2^3) (bl^3) (h^2)+12 (6^(1/3)) al aw bb0 (bb2^3) (bl^2) bw (h^2)+4 (6^(1/3)) ab0 aw (bb2^3) (bl^3) bw (h^2)+6 (6^(1/3)) (al^2) bb0 (bb2^3) bl (bw^2) (h^2)+6 (6^(1/3)) (aw^2) (bb2^2) (bl^2) (bw^2) (h^2)+6 (6^(1/3)) ab0 al (bb2^3) (bl^2) (bw^2) (h^2)+8 (6^(1/3)) al aw (bb2^2) bl (bw^3) (h^2)+(6^(1/3)) (al^2) (bb2^2) (bw^4) (h^2)+6 (6^(1/3)) al (aw^2) bb0 (bb2^3) (bl^2) (h^3)+2 (6^(1/3)) ab0 (aw^2) (bb2^3) (bl^3) (h^3)+12 (6^(1/3)) (al^2) aw bb0 (bb2^3) bl bw (h^3)+4 (6^(1/3)) (aw^3) (bb2^2) (bl^2) bw (h^3)+12 (6^(1/3)) ab0 al aw (bb2^3) (bl^2) bw (h^3)+2 (6^(1/3)) (al^3) bb0 (bb2^3) (bw^2) (h^3)+12 (6^(1/3)) al (aw^2) (bb2^2) bl (bw^2) (h^3)+6 (6^(1/3)) ab0 (al^2) (bb2^3) bl (bw^2) (h^3)+4 (6^(1/3)) (al^2) aw (bb2^2) (bw^3) (h^3)+6 (6^(1/3)) (al^2) (aw^2) bb0 (bb2^3) bl (h^4)+(6^(1/3)) (aw^4) (bb2^2) (bl^2) (h^4)+6 (6^(1/3)) ab0 al (aw^2) (bb2^3) (bl^2) (h^4)+4 (6^(1/3)) (al^3) aw bb0 (bb2^3) bw (h^4)+8 (6^(1/3)) al (aw^3) (bb2^2) bl bw (h^4)+12 (6^(1/3)) ab0 (al^2) aw (bb2^3) bl bw (h^4)+6 (6^(1/3)) (al^2) (aw^2) (bb2^2) (bw^2) (h^4)+2 (6^(1/3)) ab0 (al^3) (bb2^3) (bw^2) (h^4)+2 (6^(1/3)) (al^3) (aw^2) bb0 (bb2^3) (h^5)+2 (6^(1/3)) al (aw^4) (bb2^2) bl (h^5)+6 (6^(1/3)) ab0 (al^2) (aw^2) (bb2^3) bl (h^5)+4 (6^(1/3)) (al^2) (aw^3) (bb2^2) bw (h^5)+4 (6^(1/3)) ab0 (al^3) aw (bb2^3) bw (h^5)+(6^(1/3)) (al^2) (aw^4) (bb2^2) (h^6)+2 (6^(1/3)) ab0 (al^3) (aw^2) (bb2^3) (h^6)-2 (6^(1/3)) (bb2^3) (bl^3) (bw^2) xx-4 (6^(1/3)) aw (bb2^3) (bl^3) bw h xx-6 (6^(1/3)) al (bb2^3) (bl^2) (bw^2) h xx-2 (6^(1/3)) (aw^2) (bb2^3) (bl^3) (h^2) xx-12 (6^(1/3)) al aw (bb2^3) (bl^2) bw (h^2) xx-6 (6^(1/3)) (al^2) (bb2^3) bl (bw^2) (h^2) xx-6 (6^(1/3)) al (aw^2) (bb2^3) (bl^2) (h^3) xx-12 (6^(1/3)) (al^2) aw (bb2^3) bl bw (h^3) xx-2 (6^(1/3)) (al^3) (bb2^3) (bw^2) (h^3) xx-6 (6^(1/3)) (al^2) (aw^2) (bb2^3) bl (h^4) xx-4 (6^(1/3)) (al^3) aw (bb2^3) bw (h^4) xx-2 (6^(1/3)) (al^3) (aw^2) (bb2^3) (h^5) xx-(((9 (bb2^4) (bl^4) (bw^4) yy+36 aw (bb2^4) (bl^4) (bw^3) h yy+36 al (bb2^4) (bl^3) (bw^4) h yy+54 (aw^2) (bb2^4) (bl^4) (bw^2) (h^2) yy+144 al aw (bb2^4) (bl^3) (bw^3) (h^2) yy+54 (al^2) (bb2^4) (bl^2) (bw^4) (h^2) yy+36 (aw^3) (bb2^4) (bl^4) bw (h^3) yy+216 al (aw^2) (bb2^4) (bl^3) (bw^2) (h^3) yy+216 (al^2) aw (bb2^4) (bl^2) (bw^3) (h^3) yy+36 (al^3) (bb2^4) bl (bw^4) (h^3) yy+9 (aw^4) (bb2^4) (bl^4) (h^4) yy+144 al (aw^3) (bb2^4) (bl^3) bw (h^4) yy+324 (al^2) (aw^2) (bb2^4) (bl^2) (bw^2) (h^4) yy+144 (al^3) aw (bb2^4) bl (bw^3) (h^4) yy+9 (al^4) (bb2^4) (bw^4) (h^4) yy+36 al (aw^4) (bb2^4) (bl^3) (h^5) yy+216 (al^2) (aw^3) (bb2^4) (bl^2) bw (h^5) yy+216 (al^3) (aw^2) (bb2^4) bl (bw^2) (h^5) yy+36 (al^4) aw (bb2^4) (bw^3) (h^5) yy+54 (al^2) (aw^4) (bb2^4) (bl^2) (h^6) yy+144 (al^3) (aw^3) (bb2^4) bl bw (h^6) yy+54 (al^4) (aw^2) (bb2^4) (bw^2) (h^6) yy+36 (al^3) (aw^4) (bb2^4) bl (h^7) yy+36 (al^4) (aw^3) (bb2^4) bw (h^7) yy+9 (al^4) (aw^4) (bb2^4) (h^8) yy+(Sqrt[3]) √((bb2^6) (((bl+al h))^6) (((bw+aw h))^6) (2 ((((bw^2)+2 ab0 bb2 bl h+2 aw bw h+(aw^2) (h^2)+2 ab0 al bb2 (h^2)+2 bb0 bb2 (bl+al h)-2 bb2 bl xx-2 al bb2 h xx))^3)+27 (bb2^2) (((bl+al h))^2) (((bw+aw h))^2) (yy^2)))))^(2/3)))/((6^(2/3)) (bb2^2) (((bl+al h))^2) (((9 (bb2^4) (bl^4) (bw^4) yy+36 aw (bb2^4) (bl^4) (bw^3) h yy+36 al (bb2^4) (bl^3) (bw^4) h yy+54 (aw^2) (bb2^4) (bl^4) (bw^2) (h^2) yy+144 al aw (bb2^4) (bl^3) (bw^3) (h^2) yy+54 (al^2) (bb2^4) (bl^2) (bw^4) (h^2) yy+36 (aw^3) (bb2^4) (bl^4) bw (h^3) yy+216 al (aw^2) (bb2^4) (bl^3) (bw^2) (h^3) yy+216 (al^2) aw (bb2^4) (bl^2) (bw^3) (h^3) yy+36 (al^3) (bb2^4) bl (bw^4) (h^3) yy+9 (aw^4) (bb2^4) (bl^4) (h^4) yy+144 al (aw^3) (bb2^4) (bl^3) bw (h^4) yy+324 (al^2) (aw^2) (bb2^4) (bl^2) (bw^2) (h^4) yy+144 (al^3) aw (bb2^4) bl (bw^3) (h^4) yy+9 (al^4) (bb2^4) (bw^4) (h^4) yy+36 al (aw^4) (bb2^4) (bl^3) (h^5) yy+216 (al^2) (aw^3) (bb2^4) (bl^2) bw (h^5) yy+216 (al^3) (aw^2) (bb2^4) bl (bw^2) (h^5) yy+36 (al^4) aw (bb2^4) (bw^3) (h^5) yy+54 (al^2) (aw^4) (bb2^4) (bl^2) (h^6) yy+144 (al^3) (aw^3) (bb2^4) bl bw (h^6) yy+54 (al^4) (aw^2) (bb2^4) (bw^2) (h^6) yy+36 (al^3) (aw^4) (bb2^4) bl (h^7) yy+36 (al^4) (aw^3) (bb2^4) bw (h^7) yy+9 (al^4) (aw^4) (bb2^4) (h^8) yy+(Sqrt[3]) √((bb2^6) (((bl+al h))^6) (((bw+aw h))^6) (2 ((((bw^2)+2 ab0 bb2 bl h+2 aw bw h+(aw^2) (h^2)+2 ab0 al bb2 (h^2)+2 bb0 bb2 (bl+al h)-2 bb2 bl xx-2 al bb2 h xx))^3)+27 (bb2^2) (((bl+al h))^2) (((bw+aw h))^2) (yy^2)))))^(1/3)))];

Try this out with nominal values :

Clear[al,bl,aw,bw,ab0,bb0,bb2,alpha];al = 5.4;  bl = 10.7;aw = 5.8;  bw =1.9;ab0 = -2.4 ; bb0 = 7.9 ;bb2 = -1.46;alpha=Pi/4;

Plot[ym4[h],{h,0,10}]

Define the maximum value of x for the plot :

Plot[xm[h],{h,0,10}]

Check the top view of the brink curve shape and its intersection with the ground :

hh=5;Plot[{bx0[y,bz[y,hh],hh],bx0[y,.8*bz[y,hh],hh],bx0[y,.6*bz[y,hh],hh],bx0[y,.4*bz[y,hh],hh],bx0[y,.2*bz[y,hh],hh],bx0[y,0,hh]},{y,0,ym4[hh]}]

Check that ybrink indeed yields the y - coord at the brink just above the point at height zzz emanating from the brink at y :

hh=5;zzz=1;Plot[ybrink[bx0[y,zzz,hh],by0[y,zzz,hh],hh]-y,{y,0,10},PlotRange->All]

zslip[x_,y_,h_]:=bz[ybrink[x,y,h],h]*(x-bx0[ybrink[x,y,h],0,h])/(bx[ybrink[x,y,h],h]-bx0[ybrink[x,y,h],0,h]); (* height of the (x,y) point on the slip face *);

Define the whole shape function and restrict attention to z ≥ 0:

ztotal[x_,y_,h_]:=If[x<bx[y,h],Max[zslip[x,y,h],0],Max[zz[x,y,h],0]] (* dune height *);

Make the toe concave using a gaussian amplitude:

ztotalexp[x_,y_,h_]:=If[x>0,Exp[-((x-Max[0,bx[0,hh]])/l[hh])^2],1]*If[x<bx[y,h],Max[zslip[x,y,h],0],Max[zz[x,y,h],0]] (* dune height *);

hh=10;Plot[{ztotal[x,0,hh],ztotalexp[x,0,hh]},{x,xm[hh],l[hh]}]

Compare the two dune shapes at top precision; remove MaxRecursion -> 5 to speed up:

hh=10;Plot3D[ztotal[x,Abs[y],hh],{x,xm[hh],l[hh]},{y,-w[hh],+w[hh]},MeshFunctions->{#3&},Mesh->10,ColorFunction->Function[{x,y,z},Hue[z]],BoxRatios->{xm[hh],w[hh],hh},MaxRecursion->5]

hh=10;Plot3D[ztotalexp[x,Abs[y],hh],{x,xm[hh],l[hh]},{y,-w[hh],+w[hh]},MeshFunctions->{#3&},Mesh->10,ColorFunction->Function[{x,y,z},Hue[z]],BoxRatios->{xm[hh],w[hh],hh},MaxRecursion->5]