Here's my incomplete attempt, hopefully it will be helpful!
By symmetry, we know that the surface of the droplet must be rotationally symmetric about some axis perpendicular to the horizontal surface, so we can describe it as a solid of revolution about this axis. We can assume that the droplet only makes contact with the horizontal surface up to some maximum value $r_0$, and has a height $h_0$ at its centre. It is best to describe the surface using the radius as a function $r(h)$ of height, as describing the height as a function of radius $h(r)$ would not allow the description of droplets which overhang $r_0$. This is shown in the figure below. I've made the droplet wobbly because we don't yet know what shape it will take.

The shape which surface will take is the one which minimises the total energy of the configuration, while keeping the volume of the drop constant. There are two contributions to the energy: the gravitiational potential energy, which depends only on the height of the centre of mass of the droplet; and the surface areas of the various interfaces, each multiplied by the constants $\gamma$. Let $S_{\mathrm{lg}}$ be the area of the curved surface of the drop, $S_{\mathrm{ls}}$ be its contact area with the solid plane, $S_{\mathrm{sg}}$ be the area of the solid surface in contact with the air, $V$ be the total volume of the drop, $\rho$ its density (assumed constant), and finally $\bar{h}$ the height of its centre of mass. Then the energy is
$$
E =\gamma_{\mathrm{lg}}S_{\mathrm{lg}} + \gamma_{\mathrm{ls}}S_{\mathrm{ls}} + \gamma_{\mathrm{sg}}S_{\mathrm{sg}} +g V\rho\bar{h}.
$$
The area of the entire plane surface is constant, so $S_{\mathrm{ls}} + S_{\mathrm{sg}} = \mathrm{const}$, so, to within a constant offset which doesn't matter, the energy is
$$
E =\gamma_{\mathrm{lg}}S_{\mathrm{lg}} + (\gamma_{\mathrm{ls}} - \gamma_{\mathrm{sg}})S_{\mathrm{ls}} + g V\rho\bar{h}.
$$
This answers the first of your questions: yes, in general energy is stored in the solid-gas interface, but it can be accounted for simply by changing the liquid-solid surface tension constant from $\gamma_{\mathrm{ls}}$ to the 'effective' surface tension constant $\gamma_{\mathrm{ls}} - \gamma_{\mathrm{sg}}$.
We now need to express the various terms in our expression for the energy as functionals of $r(h)$. It is not too difficult to see that
\begin{align}
\rho V\bar{h}[r(h)] & = \int_0^{h_0}\rho \pi r^2 h \,\text{d}h,\\
V[r(h)] & = \int_0^{h_0} \pi r^2 \,\text{d}h,\\
S_{\mathrm{lg}}[r(h)] &= \int_0^{h_0} 2\pi r\sqrt{1+(r')^2} \,\text{d}h,\\
S_{\mathrm{ls}}[r(h)] &= \pi r(0)^2 = \int_0^{h_0} \pi r^2 \delta(h) \,\text{d}h.
\end{align}
In the last line I have used the Dirac delta 'function' to express $S_{\mathrm{ls}}$ as an integral - I'm not sure if this is a sensible thing to do or not.
To minimise the energy functional while keeping the volume constant, we use Lagrange multipliers and so must minimise the functional
\begin{align}
I[r] &= E[r] +\lambda V[r]\\
&=\pi\int_0^{h_0} \underbrace{\left[2\gamma_{\mathrm{lg}}r\sqrt{1+(r')^2} +(\gamma_{\mathrm{ls}} - \gamma_{\mathrm{sg}})r^2 \delta(h) + g\rho r^2h +\lambda r^2\right]}_{\mathcal{L}(h, r, r')}\,\text{d}h.
\end{align}
Using the Euler-Lagrange equation
$$
\frac{\partial\mathcal{L}}{\partial r} = \frac{\text{d}}{\text{d} h}\frac{\partial\mathcal{L}}{\partial r'}
$$
gives us (after some simplification)
$$
g\rho r h +\gamma_{\text{lg}}\sqrt{1+(r')^2} +(\gamma_{\text{ls}}-\gamma_{\text{sg}})r \delta +\lambda r = \frac{\gamma_{\text{lg}}}{(1+(r')^2)^{3/2}}((r')^2+(r')^4 +r''r).
$$
In theory, solving this horrible looking differential equation should give the correct solution! Since there is a delta function, I'm not even sure if it's well-posed. I would be interested to hear if anyone knows how to proceed from here.
Note that setting $g =0$ in the above equation gives solutions of the form $r = \sqrt{\frac{1}{\lambda^2\gamma_{\text{lg}}^2}- (h-b)^2}$ for $h>0$, where $(h_0-b)^2=1/\lambda\gamma_{\text{lg}}$. This shows that when gravity is neglected the shape of the droplet is a sphere.
Once you have solved the above equation, you can use the constraint on the volume to determine the constant $\lambda$. The functional form of $r(h)$ would also allow you to find the wetting angle, using $\tan\alpha =r'(0) $.
Finally, we note that the endpoint $h_0$ is variable. This means that, in addition to the Euler-Lagrange equations, we require that
$$
\left.\mathcal{L} - r'\frac{\partial\mathcal{L}}{\partial r'} \right|_{h=h_0} = 0
$$
(see e.g. here for a derivation). This condition says that
$$
\left.2\gamma_{\mathrm{lg}}r\sqrt{1+(r')^2} +(\gamma_{\mathrm{ls}} - \gamma_{\mathrm{sg}})r^2 \delta(h) + g\rho r^2h +\lambda r^2 - r'\cdot 2\gamma_{\text{lg}}r'r\sqrt{1+(r')^2} \right|_{h=h_0} = 0
$$
which is true when $r(h_0) = 0$. This condition is therefore equivalent to the height of the droplet being $h_0$.
ground
in this context is misleading. – Alex Trounev Apr 07 '20 at 10:50