Instantaneity of Stokes Flow Background Info:
A steady-state incompressible newtonian Stokes fluid occupying a domain $\Omega \in \mathbb{R}^2$ has flow velocity $\textbf{u}$ and pressure $p$ given by the following system of PDEs:
\begin{align*}
    -\nu \Delta \textbf{u} + \frac{1}{\rho} \nabla p &= \textbf{f} \text{ in } \Omega, \\ 
    \nabla \cdot \textbf{u} &= 0 \text{ in } \Omega. \label{eq:4}
\end{align*}
In order for this system to have a unique solution, values for $\textbf{u}$ (or its derivatives) must be prescribed at the boundaries, and some condition on $p$ must be enforced in order to overcome certain uncertainty (since $p$ appears through its gradient and is thus only determined up through an arbitrary constant). I am interested in solving this system of PDEs with (time-dependent) moving boundary conditions numerically (using a Finite Element Method). I am unsure if I should use the steady-state versus unsteady equations that describe the motion of Stokes fluid.

More relevant information directly pertaining to my question:
The wikipedia article on Stokes flow states that Stokes fluids are instantaneous, meaning there is no dependence on time in the system above (other than possibly through boundary conditions). However, the article also states that "occasionally one might consider the unsteady Stokes equations, in which the term $ \rho {\frac{\partial \mathbf {u}}{\partial t}}$ is added to the left hand side of the momentum balance equation". 

Specific questions:
In this case, wouldn't adding the time-dependent term $\rho \frac{\partial \textbf{u}}{\partial t}$ make the Stokes fluid non-instantaneous? Isn't the point of a "steady-state" fluid to be time-independent anyways? Why is it even mentioned that Stokes flow is instantaneous then? What exactly does adding the time-dependent term do if "given the boundary conditions of a Stokes flow, the flow can be found without knowledge of the flow at any other time." If I want to simulate Stokes flow in a box, where one of the boundaries is moving, do I need to include the $\rho \frac{\partial \textbf{u}}{\partial t}$ term?.
 A: If the question is "Do I need to include the $\rho \frac{\partial \textbf{u}}{\partial t}$ term to simulate Stokes flow", the answer is no because (as you point out) that term is neglected in the definition of the assumption one makes for obtaining Stokes flow.
The Wikipedia article you cite is using a poor choice of words to state that sometimes researchers want to model non-instantaneous effects in low-Reynolds hydrodynamics while retaining a lot of the "infrastructure" associated with the Stokes equations, namely neglecting the convective acceleration. Some people call that "time-dependent Stokes equations" or "transient Stokes equations" because of the similarity, but it leads to misunderstandings regarding the physical consequences of the assumptions in "mainline" Stokes flow, as you've pointed out. 
A: "If I want to simulate Stokes flow in a box, where one of the boundaries is moving, do I need to include the $\rho \frac{\partial \mathbf{u} }{\partial t}$ term?."
One possible way to answer this question is to simply examine the dimensionless form of the Navier-Stokes equations.  You want to essentially understand whether the intertial acceleration term is negligible.  For low Reynolds number flows the inertial convective $\mathbf{u} \cdot \nabla \mathbf{u}$ term drops (assuming $Re \ll 1$) and the non-dimensional Navier-Stokes equations reduce to,
$$
\left( \frac{L^{2}}{\nu T} \right) \frac{ \partial \bar{\mathbf{u} }}{\partial \bar{t} } = -\bar{\nabla}\left( \bar{P}\right)+ \bar{\nabla}^{2}\bar{\mathbf{u}} + \frac{Re}{Fr}\frac{\mathbf{f}}{|\mathbf{f}|}
$$
where $\bar{\mathbf{x}} = \mathbf{x}/L$, $\bar{\mathbf{u}} = \mathbf{u}/U$, $\bar{t} = t/T$, $\bar{\nabla} = \frac{\partial}{\partial \bar{\mathbf{x}}}$, $\bar{\nabla}^{2} = \frac{\partial}{\partial \bar{x}_{i} \bar{x}_{i}}$, $Re = UL/\nu$, and $Fr = U^{2}/|\mathbf{b}|L$.
Here, $L$ is a characteristic length related to the size of the boundaries, $U$ is a characteristic velocity determined by the mechanism driving the flow, and $T$ is a characteristic time either imposed by external forcing or simply defined by $L/U$.  In the latter case you recover the Reynolds number in from of the inertial acceleration piece.
For your moving boundary, can you calculate $L^{2}/\nu T$?  This quantity expresses the magnitude of inertial acceleration forces relative to viscous forces.  In this expression, $L$ would be the boundary length, $\nu$ is a fluid property (e.g. you are simulating water at STP conditions use an appropriate density and dynamic viscosity), and $T$ is the characteristic time associated with your boundary forcing mechanism.  This will at least tell you how this term scales and allow you to judge whether you need to account for it.
As a side note, the solution of the Stokes equations gives a time-independent flow field solution that propagates instantaneously in space.  This is a fine assumption so long as there is reason to drop the relevant terms shown above.  This approach is used extensively to simulated the hydrodynamics of small active particles (the field of Active Matter) and suspensions of particles.  Usually when one wants to understand the motions of objects immersed in a Stokesian fluid, the Mobility Problem is solved which yields rigid body motions (linear and angular velocities of particles) on specification of applied external forces.  The Stokes equations are solved to find particle velocities and these velocities are used to time step (e.g. using the Euler method) the physical location of the particles and map their trajectories.
