Footings Under Biaxial Bending
When designing concrete spread footings, finding its strength part is rather simple, using basic principles of concrete design. Finding the loads to design the footing however can be much trickier, particularly in the case where there is bending in both the X and Y direction simultaneously. This article will go into detail about how ClearCalcs performs these calculations  and how you can do the same by hand if you're feeling adventurous.
First, let us establish the four main dimensions which the rest of our calculations will depend on. Our footing plan dimensions are defined by the width $B$, parallel to the X axis and the length $L$, parallel to the Y axis.
The relationships between moment and axial load is best defined by the eccentricities of the load, defined as such:
Note that for ease of calculation, we take the absolute values of moments, which means that the resultant force will always be in the first quadrant. Because of symmetry, the values we get below will all be the same regardless of which quadrant this occurs in.
Soil Bearing Calculations
These equations are all provided by John Bellos and Nikolaos Bakas in their paper available here, and were independently verified by our team. They are based on first principles of statics and satisfy equilibrium of forces and moments for all cases. The symbol convention we use differs slightly from the paper, but the equations remain the same.
Eccentricity Zones
The most difficult part of finding the stress distribution is that since soil cannot resist uplift by tension, at high eccentricities some parts of the footing may lift off the ground and not carry any load. The shape of the area still carrying load will vary based on the eccentricity of the loads, and will correspond to one of the five shapes shown below.
Zone  Shape of Loaded Area  Eccentricity Limits 

1 

$$\dfrac{e_x}{L} + \dfrac{e_y}{B} \leq \dfrac{1}{6}$$ 
2  $$\dfrac{e_x}{L} + \dfrac{e_y}{B} \gt \dfrac{1}{6}$$ and $$e_x \lt \dfrac{L}{2}  \dfrac{L}{6} \left(2\sqrt{112\dfrac{e_y^2}{B^2}}\right)\left(1+2\dfrac{e_y}{B}+\sqrt{112\dfrac{e_y^2}{B^2}}\right) \Big/\left(1 + 4\dfrac{e_y}{B}\right)$$ and $$e_y \lt \dfrac{B}{2}  \dfrac{B}{6} \left(2\sqrt{112\dfrac{e_x^2}{L^2}}\right)\left(1+2\dfrac{e_x}{L}+\sqrt{112\dfrac{e_x^2}{L^2}}\right) \Big/\left(1 + 4\dfrac{e_x}{L}\right)$$  
3  $$e_y \geq \dfrac{B}{2}  \dfrac{B}{6} \left(2\sqrt{112\dfrac{e_x^2}{L^2}}\right)\left(1+2\dfrac{e_x}{L}+\sqrt{112\dfrac{e_x^2}{L^2}}\right) \Big/\left(1 + 4\dfrac{e_x}{L}\right)$$ and $e_x\leq \dfrac{L}{4}$  
4  $$e_x \geq \dfrac{L}{2}  \dfrac{L}{6} \left(2\sqrt{112\dfrac{e_y^2}{B^2}}\right)\left(1+2\dfrac{e_y}{B}+\sqrt{112\dfrac{e_y^2}{B^2}}\right) \Big/\left(1 +4\dfrac{e_y}{B}\right)$$ and $e_y\leq \dfrac{B}{4}$  
5  $e_x\gt \dfrac{L}{4}$ and $e_y\gt \dfrac{B}{4}$ 
These equations may be easy for a computer to digest, but they are a handful to look at, particularly for zones 2 to 4. Based on the eccentricity limits shown above, we can map five different zones on the footing where, if the eccentricity of the load falls within this area, the shape of the loaded area will correspond to the shapes shown above. This map of the eccentricity zones makes it much easier to see where the eccentricity falls. We generated this map using the Desmos online graphical calculator  see it here!
Zones 3, 4 and 5 extend quite far from the center of the footing  in particular, the overturning factor of safety will always be very low if the eccentricity is in zone 5, and the design will be far from efficient.
Bearing Pressure Profile
Once we've identified the eccentricity zone, we know the shape that our bearing area will take. At this point, we now need to define it. In their paper, Bellos & Bakas define the bearing area with 3 points  the maximum bearing pressure ($q_{max}$), and the two intercepts of the zero pressure line (the neutral axis of the footing) with the axes formed by the two upperright edges of the footing ($x_n$ and $y_n$). The intercepts can be outside the footing (low eccentricity) or inside (high eccentricity).
This method proves quite convenient later when it comes time to find the footing loads from the bearing pressure, as will be shown later. The three parameters, $q_{max}$, $x_n$ and $y_n$ are defined in the sections below. These were derived by Bellos & Bakas using principles of equilibrium and statics.
Once we've got our bearing profile properties, we can compare the maximum bearing pressure with our allowable bearing pressure. This is what will govern the plan dimensions of the footing. Once these dimensions are determined, we move on to determine the shear and moment loads on the critical footing section, which we will go over in the next section.
Zone 1
This is the zone with the lowest eccentricity, where the entire footing is under compression. The intercepts are given by:
and
The maximum bearing pressure:
Zone 2
This is the troublesome zone, where there is no closedform solution  this case requires solving a system of two 6th order polynomials  so we must use a numerical solving method (ClearCalcs makes use of the multivariate Newton's method). First we solve for the two intercepts using the system of equation below:
and
Using the Newton method, this converges extremely quickly  3 iterations is almost always enough to get less than 1% error. Once this is solved, we can find $q_{max}$:
Zone 3
This is the case where eccentricity in the $y$ direction is much larger than in the $x$ direction. Finding the intercepts:
and
We can then find the maximum bearing pressure:
Zone 4
This is the same as Zone 3, except that the $x$ and $y$ eccentricities are reversed. Finding the intercepts:
and
We can find the maximum bearing pressure:
Zone 5
This is perhaps the easiest zone to calculate, and these equations can easily be derived graphically by using properties of tetrahedrons. Finding the intercepts:
and
The maximum bearing pressure is also found:
Footing Loads
Finding the loads acting on the concrete footing is the next step in our design  this is used to determine the thickness of the footing and the reinforcement required. The general basis of the procedure is that we must find the volume under the bearing pressure profile to find the total shear load, and its centroid to find the bending moment.
In general, most software packages make use of double or triple integration which is performed numerically. This however is relatively slow and more importantly, is very "blackboxy", which we try to avoid as much as we can. Thus, we derived our own method of finding shear and moment loads on the footing, which we will explain below.
Using Tetrahedrons to Represent the Bearing Pressure Profile
If we look at the bearing profile below, we can clearly see that what we have is a tetrahedron (tet), of height $q_{max}$ and a base of dimensions $x_n \times y_n$.
This is a special case of a tet with three right angles at the same vertex  this is called a trirectangular tetrahedron. Its volume and centroid are both easy to find:
The problem arises when our tetrahedron extends past the width of the footing or past the critical section, which we define as a distance $a$ from the edge of the footing. All of a sudden, our nicely defined tet gets chopped up and these easy properties don't apply so well. In this case, we need to split our tet into different components, as shown below:
We can now see that we have up to four different individual tets which are also trirectangular, labelled from A to D. Tet A contains both B, C and D, while Tet B and Tet C both contain Tet D. This becomes important when we find our footing loads later. Their properties depend on the soil bearing pressure profile, the footing dimensions, and the critical section location, which is normally different for shear and moment.
For shear:
And for moment:
Where $\ell$ is the column length.
Tet A
This is the simplest one to find  it is the ideal tet we defined above. The only difference is that we now define our centroid by its distance from the critical section, whereas it was previously defined by its distance from the footing edge.
Tet B
This tet is associated with a lower value of $e_y$. It shows up when:
The width is easy to find, but to find the length and height, we use similar triangles. We obtain our volume:
To find our centroid, we know that it is at one quarter of the length of the tet away from the footing edge. Thus the centroid about the critical section:
Tet C
This tet is associated with a lower value of $e_x$. It shows up when:
In this case, the length is easy to find, but to find the width and height, we use similar triangles. We obtain our volume:
We find the centroid about the critical section, but in this case, it is below the critical section so we add a negative sign:
Tet D
This tet only shows up in cases of low eccentricities, when the length of Tet B is greater than the distance to the critical section:
This is by far the trickiest one to find, stuck inside both Tet B and Tet C. But this also means we can use some of their properties to make our life easier. For the volume:
Finding the centroid, which is also below the critical section:
Finding the actual footing loads
Now that we've got the properties of all our tets, we can go back to finding the shear and moment loads. This is actually the easy part  shear is just the combination of volumes and moment is the combination of volumes multiplied by their moment arms! Bringing back our diagram:
We can come up with a general equation to find the shear:
And for moment:
Depending on the values of $x_n$ and $y_n$, some of the tets may disappear in which case their associated terms are removed from the previous equations. There are five different possibilities:

$x_n\left(1\frac{B}{y_n}\right) > a$ and $y_n > B$ 

$x_n\left(1\frac{B}{y_n}\right) \leq a < x_n$ and $y_n > B$ 
$x_n \gt a$ and $y_n \leq B$  

$x_n \leq a$ and $y_n \gt B$ 
$x_n \leq a$ and $y_n \leq B$ 
And with this, we've successfully managed to determine both our bearing pressures and the loads on our footing.