In practice, we approximate the operator \({\rm {std}}_{\omega _f,\omega _g}\equiv ({\rm {var}}_{\omega _f,\omega _g})^{0.5}\)  by the two-dimensional spherical cubature integration [1]} with 4 quadrature points (although extension to more quadrature points is straightforward):
\(m = \mathbb {E}_{\omega _f, \omega _g}[\mathbb {P} & (S(x, \omega | \omega _f, \omega _g) > \delta )] = \frac{1}{4} \sum _{i=1}^{4} \mathbb {P}(S(x, \omega | u_i )),  \\{\rm {var}}_{\omega _f, \omega _g}[\mathbb {P} & ((S(x, \omega | \omega _f, \omega _g) > \delta )] = \frac{1}{4} \sum _{i=1}^{4} (\mathbb {P}(S(x, \omega | u_{i} )) - m)^2,  \\u_1 = & \lbrace f(x) = \mu _f(x) + \sqrt{2 \; {\rm {cov}}_f(x)}, \;g(x) = 0\rbrace  \nonumber \\u_2 = & \lbrace f(x) = 0, \; g(x) = \mu _g(x) + \sqrt{2 \; {\rm {cov}}_g(x)}\rbrace  \nonumber \\u_3 = & \lbrace f(x) =-\mu _f(x) - \sqrt{2 \; {\rm {cov}}_f(x)}, \;g(x) = 0\rbrace  \nonumber \\u_4 = & \lbrace f(x) = 0, \; g(x) = - \mu _g(x) - \sqrt{2 \; {\rm {cov}}_g(x)}\rbrace  ,\) 
