Fix $a,b \ge 1$. Starting with a deck of $a+b$ distinct cards in order, consider cutting the deck into the first $a$ cards (collectively $A$) and the last $b$ cards (collectively $B$), and then riffling the two decks together (ie, merging them while preserving the order of each deck). There are ${a+b \choose a,b}$ possible outcomes. Each outcome has a number of inversions. When selecting an outcome uniformly at random, what is the distribution of the number of inversions?
Here are a couple equivalent formulations of this:
Form a grid graph with $a$ vertical steps and $b$ horizontal steps, with all edges having unit length and oriented down or right. The riffles are naturally in bijection with the paths from the top left to the bottom right: whenever the $i$-th step goes down, the $i$-th card in the riffle comes from $A$, and whenever the $i$-th step goes right, the $i$-th card comes from $B$. The number of inversions appears as the area below this path. So this question is, equivalently, what is the distribution of the area below a uniformly random path in the grid?
Select, uniformly at random, non-negative integers $\lambda_0,\ldots,\lambda_b$ subject to $\lambda_0+\cdots+\lambda_b = a$, and output $\sum_{i=0}^b i\cdot \lambda_i$. This is distributed the same as above: each choice of $\lambda_0,\ldots,\lambda_b$ corresponds bijectively to a choice of path in the grid interpretation where $\lambda_i$ counts the number of vertical steps following the $i$-th horizontal step. The sum counts, row-by-row, the area below the path. One can of course swap the roles of $a$ and $b$ as well.
Form the polynomial $$\sum_{d=0}^{ab} c_{a,b,d}\cdot q^d = \binom{a+b}{a}_q$$ in the indeterminate $q$. (See Gaussian binomial coefficient.) Sample the value $d$ with probability $c_{a,b,d}/\binom{a+b}{a}$. (This is a probability distribution.)
In any case, let's call the random variable $V_{a,b}$. It is, of course, integer-valued. Every integer in $[0,ab]$ is a possible outcome. By the 180° rotational symmetry of the grid, the distribution is symmetric about $ab/2$. In particular the mean is $ab/2$. I've computed the variance as $\frac{1}{12}ab(a+b+1)$. That's essentially the limit of what I can prove though. Empirically, the distribution is unimodal, and no outcome is particularly likely (maximum probability of any one outcome bounded by $O(1/{\sqrt{ab(a+b)}})$). In general, it seems to loosely resemble a normal distribution with the same first two moments, even when $a$ and/or $b$ is small. I'd like to understand to what extent that is true.
In case it helps, here are two more tidbits:
$V_{a,b}$ satisfies a recurrence in the sense that the distribution of $V_{a,b}$ coincides with the distribution where, with probability $\frac{a}{a+b}$, one samples $V_{a-1,b}$, while with probability $\frac{b}{a+b}$, one samples $V_{a,b-1}$ and adds $a$.
It seems nicer to work with a related quantity, $\delta_{a,b} = ab - 2V_{a,b}$. In the deck-shuffling formulation, it counts $-1$ for each inversion, and $+1$ for each non-inversion. In the grid formulation, it counts the area above the path minus the area below the path. Similarly there is an interpretation in the integer-partition formulation. $\delta_{a,b}$ is integer-valued, and every integer in $[-ab,ab]$ with the same parity as $ab$ is a possible outcome. It has a similar recurrence to $V_{a,b}$ in that $\delta_{a,b}$ coincides with the distribution where, with probability $\frac{a}{a+b}$, one samples $b+\delta_{a-1,b}$, and with probability $\frac{b}{a+b}$ one samples $-a+\delta_{a,b-1}$. $\delta_{a,b}$ is symmetrically distributed about 0 (so all odd-order moments are zero), and its variance is $\frac{1}{3}ab(a+b+1)$. Its fourth moment is $$\frac{1}{15}ab(a+b+1)\cdot[5(a^2b+b^2a) - 2(a^2+b^2) + 3ab - 2(a+b)].$$ I computed the second and fourth moments using the recurrence for $\delta_{a,b}$.