Here what I did was took the distances between the plane and multiplied them getting option A. I don't know the correct answer, but my method seems flawed, as distance between any two parallel planes will give me the height, but the area of the base can't be defined as the multiplication of two planes distance, as the base will be a parallelogram and not a rectangle. How to approach this question?
One way to go is to find a set of coterminal edges of the paralellepiped and then compute the scalar triple product (a.k.a. the determinant) of the vectors that represent these edges, per Nick Pavlov’s comment.
Another approach is to convert to a coordinate system in which the parallelepiped becomes a rectangular box. The volume of the original paralellepiped is the volume of this box times the Jacobian determinant of the coordinate transformation. This involves inverting a $4\times4$ matrix, but one of the rows is very simple, so this computation isn’t quite as bad as it sounds.
One way to compute this inverse is via row-reduction: We’ll map the first plane of each parallel pair to the coordinate planes, so form the augmented matrix $$\left[\begin{array}{rrrr|rr}1&3&-1&-2 & 1&0&0&0 \\ 2&-1&1&0 & 0&1&0&0 \\ 1&1&1&3 & 0&0&1&0 \\ 0&0&0&1 & 0&0&0&1\end{array}\right]$$ in which the first three rows are the homogeneous vectors that represent the planes. Row-reducing this matrix results in $$\left[\begin{array}{rrrr|rr}1&0&0&0 & \frac14&\frac12&-\frac14&\frac54 \\ 0&1&0&0 & \frac18&-\frac14&\frac38&-\frac78 \\ 0&0&1&0 & -\frac38&-\frac14&\frac78&-\frac{27}8 \\ 0&0&0&1 & 0&0&0&1\end{array}\right].$$ The right-hand side of this matrix is the desired coordinate transformation. Now, transform the other three planes. In theory, this requires computing the matrix product $$\left[\begin{array}{r}1&3&-1&3\\2&-1&1&3\\1&1&1&-3\end{array}\right]\left[\begin{array}{r}\frac14&\frac12&-\frac14&\frac54 \\ \frac18&-\frac14&\frac38&-\frac78 \\ -\frac38&-\frac14&\frac78&-\frac{27}8 \\ 0&0&0&1\end{array}\right]$$ but by construction we know that the result will be of the form $[I_3\mid\mathbf b]$, so we really only need to compute the last column, which consists of the dot products of the rows of the left-hand matrix with the last column of the right-hand matrix. The images of the first three planes all have zero for this last coordinate, so the absolute values of these dot products are the edge lengths of the box.
All that’s left to do is to find the determinant of the coordinate transformation. That looks a bit messy, but recall that it’s the reciprocal of the determinant of the matrix that was inverted. That looks much easier to compute. The last row consists of zeros and a one, so the determinant is equal to the upper-left minor, and that has a row of ones, so it’s just the sum of three $2\times2$ determinants.
These two considerations suggest a further optimization. Since we don’t need the entire transformation matrix, there was no need to augment with the full identity matrix. Using only its last column suffices.
This seems like a lot of work, but I think it’s less work than the first method, which involves solving four systems of three linear equations each and then computing a $3\times3$ determinant.
End note: It might seem like the transformation “scale factor” is inverted, but the vectors that represent the planes are covariant.