How do you model a rolling shutter in the frequency domain
A rolling shutter introduces a displacement of the lines of an image that is proportional to the relative velocity of the camera (or target) and the velocity of the shutter. In terms of "simulating" this through a video, it translates to "borrowing" a few lines of pixels from each frame to compose one final frame. Notice here how the shearing produced by the rolling shutter, which could be seen as a displacement of $x$ in 2 directions just got translated to a simple "vertical" motion in 3 dimensions. 

So, all that we have to do is look at what exactly this operation does in the time domain and see how (if possible) we can replicate it in the frequency domain.

From your question, I suppose that you want to create one image (let's call it $I$) from the "integration" of the frames in $a$. Before we construct $I$, we have to make some assumptions about how $a$ was acquired. So, let's assume that $a$ was acquired with a "perfect" camera. One that does not suffer from rolling shutter itself and one that shoots at a very high frame rate. So the period of sampling $T_s = t_k - t_{k-1}$ of $a$ is $T_s \ll T_e$, where $T_e$ is the frame exposure time (the time it takes for the shutter to roll and create one image). Furthermore, let's assume that the "width" of the shutter, in terms of image lines, is $W_s$. 

Then, just considering the shutter, this is indeed a "white" rectangle that moves parallel to axis $y$ across all time instances of $a$ and it looks like:

$$S(x,y,t) = rect\left(\frac{x-\frac{N_x}{2}}{\frac{N_x}{2}},\frac{y-t \times T_e}{\frac{W_s}{2}}\right)$$

This is a two dimensional rectangle. It is centred in the middle of the $X$ dimension of the frame and is as wide as the frame and in the $Y$ direction it is $W_s$ lines tall, moving downwards.

The time instances that make up $I$ are:

$$I(x,y,t) = a(x,y,t) \times S(x,y,t)$$

But if you "replay" this, all that you see is a band of pixels from each frame moving downwards.

To get rid of $t$ in $I$, you need to integrate this signal with respect to $t$. In other words, you have to sum all pixel time courses to generate just one image that is independent of $t$. 

Now, this is a linear time invariant system of operations in three dimensions. To find its coefficients, you need to first transform $S_{FT} = \mathscr{F}(S(x,y,t))$.

To produce $I$ we multiplied $a$ with $S$. Multiplication in the time domain is convolution in the frequency domain. Therefore to express $I(x,y,t)$ (notice the $t$ still there):

$$I_{FT}(k,m,f) = \mathscr{F}(a(x,y,t)) \ast \mathscr{F}(S(x,y,t))$$

Here of course, notice the change of variables, $m,n,f$ represent their frequency domain $x,y,t$ counterparts.

Right, so now, we want to integrate the three dimensional complex $I_{FT}(k,m,f)$ with respect to time. We look up the properties of the Fourier Transform and realise that integration in time is equivalent to multiplication in frequency by a $\frac{1}{j \omega}$ factor plus a static DC component scaled by $\pi \delta(\omega)$. So, all we have to do now is multiply our $I_{FT}$ coefficients and offset them by these factors, in the direction of the time vector. Or in other words, in the dimension that corresponds to frames.

$$I_{FT}(k,m,f) = I_{FT}(k,m,f) \times \frac{1}{e^{2 \pi i f/N_f}} + \pi I_{FT}(0,0,0)$$

(And of course then, above, obtain the inverse to get back into "time" domain)

Of course, you can express the two as one operation and do the transformation in one step.

Hope this helps.