On a sphere of radius $R > 0$, a geodesic triangle with interior angles $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ has area $R^{2}(\theta_{1} + \theta_{2} + \theta_{3} - \pi)$. One way you might proceed, therefore, it to triangulate your geodesic polygon (whatever its actual shape) and sum the areas of its triangular pieces. For a quadrilateral with interior angles $\theta_{1}$, $\theta_{2}$, $\theta_{3}$, and $\theta_{4}$, the area is
$$
R^{2}(\theta_{1} + \theta_{2} + \theta_{3} + \theta_{4} - 2\pi).
$$
Added (in light of OP's clarifications): If $\theta \leq \pi$ is the interior angle at each vertex of the quadrilateral, then
$$
\theta = \pi - \cos^{-1}\left\lvert\tan\frac{\ell}{2R} \tan\frac{b}{2R}\right\rvert,
$$
so
$$
\text{Area}
= R^{2}(4\theta - 2\pi)
= R^{2} \left[2\pi - 4\cos^{-1}\left\lvert\tan\frac{\ell}{2R} \tan\frac{b}{2R}\right\rvert\right].
$$
(Particularly, the absolute value of the product of tangents inside the arccos does not exceed unity.)
To see this, it's convenient to work in Cartesian coordinates with the sphere centered at the origin. Denote the vertices of the quadrilateral by
$$
v_{1} = (A, B, C),\qquad
v_{2} = (A, -B, C),\qquad
v_{3} = (A, -B, -C),\qquad
v_{4} = (A, B, -C).
$$
Of course, $A^{2} + B^{2} + C^{2} = R^{2}$.

The plane through the origin, $v_{1}$, and $v_{2}$ has equation $Cx - Az = 0$, and so has unit normal vector
$$
n_{1} = \frac{(C, 0, -A)}{\sqrt{A^{2} + C^{2}}} = \frac{(C, 0, -A)}{\sqrt{R^{2} - B^{2}}}.
$$
The plane through the origin, $v_{1}$, and $v_{4}$ has equation $Bx - Ay = 0$, and so has unit normal vector
$$
n_{2} = \frac{(B, -A, 0)}{\sqrt{A^{2} + B^{2}}} = \frac{(B, -A, 0)}{\sqrt{R^{2} - C^{2}}}.
$$
The "large" angle between these planes is the "large" angle between the great circles they determine (because each normal vector $n_{i}$ is tangent to the sphere at $v_{1}$), i.e., the interior angle $\theta$ of the quadrilateral. Taking the ordinary dot product of the unit normals,
$$
\cos\theta = n_{1} \cdot n_{2}
= \frac{B}{\sqrt{R^{2} - B^{2}}}\, \frac{C}{\sqrt{R^{2} - C^{2}}}.
\tag{1}
$$
Let $2\psi_{1}$ denote the angle subtended at the center of the sphere by the side from $v_{1}$ to $v_{2}$, and let $\ell_{1} = 2R\psi_{1}$ denote the corresponding side length, so that $\psi_{1} = \ell_{1}/(2R)$. Thinking of $v_{1}$ and $v_{2}$ as vectors in space, the ordinary dot product gives
$$
\frac{v_{1} \cdot v_{2}}{R^{2}}
= \cos(2\psi_{1})
= 1 - 2\sin^{2} \psi_{1}.
$$
On the other hand, using the Cartesian components of these vectors, we have
$$
\frac{v_{1} \cdot v_{2}}{R^{2}}
= \frac{A^{2} - B^{2} + C^{2}}{A^{2} + B^{2} + C^{2}}
= 1 - 2\frac{B^{2}}{R^{2}}.
$$
Equating, it follows at once that $B/R = \pm\sin\psi_{1}$. (As a consistency check, the isoceles triangle with the origin, $v_{1}$, and $v_{2}$ as vertices has apex angle $2\psi_{1}$ and base $2B$.)
Similarly, letting $2\psi_{2}$ denote the angle subtended at the center of the sphere by the side from $v_{1}$ to $v_{4}$, and letting $\ell_{2} = 2R\psi_{2}$ denote the corresponding side length, we find $C/R = \pm\sin\psi_{2}$.
That is,
$$
B = \pm R\sin\frac{\ell_{1}}{2R},\qquad
C = \pm R\sin\frac{\ell_{2}}{2R},
\tag{2}
$$
and consequently
$$
\frac{B}{\sqrt{R^{2} - B^{2}}} = \pm\tan\frac{\ell_{1}}{2R},\qquad
\frac{C}{\sqrt{R^{2} - C^{2}}} = \pm\tan\frac{\ell_{1}}{2R}.
\tag{3}
$$
Substituting (3) into (1),
$$
\cos\theta
= \frac{B}{\sqrt{R^{2} - B^{2}}}\, \frac{C}{\sqrt{R^{2} - C^{2}}}
= \pm \tan\frac{\ell_{1}}{2R} \tan\frac{\ell_{2}}{2R}.
$$
To get the large angle, either take the arccos of the negative value, or $\pi$ minus the arccos of the positive value. (The angle/area formulas above do the latter.)