I'm going to change your notation for (3), so that my answer will be compact enough to fit on one page.
Restating your model, $$\vec x_k = \left [ \begin{array}{c|c} \mathbf 1 & \mathbf T \\ \hline \mathbf 0 & \mathbf 1 \end{array} \right ] x_{k-1} + \begin{bmatrix}\mathbf{\frac {T^2}{2}} \\ \hline \mathbf T \end{bmatrix} u_k \tag a$$
This says the same thing, but leaves you to infer that the elements are all 3x3 identity matrices multiplied by the given factor. Note that I've trimmed your $u_k$ to just three elements to eliminate redundancy.
To add accelerometer bias into the mix, just add it as a state, so that your (1) becomes $$ \vec x = \begin{bmatrix} p_x & p_y & p_z & v_x & v_y & v_z & \Delta a_x & \Delta a_y & \Delta a_z \end{bmatrix}^T \tag b$$ with $\Delta a_x \cdots$ being the x, y, and z components of the accelerometer bias.
Now augment your model for the extra states: $$\vec x_k = \left [ \begin{array}{c|c|c} \mathbf 1 & \mathbf T & \mathbf{\frac{T^2}{2}} \\ \hline \mathbf 0 & \mathbf 1 & \mathbf T \\ \hline \mathbf 0 & \mathbf 0 & \mathbf 1 \end{array} \right ] x_{k-1} + \begin{bmatrix}\mathbf{\frac {T^2}{2}} \\ \hline \mathbf T \\ \hline \mathbf 0 \end{bmatrix} u_k \tag a$$
Then choose reasonable numbers for your accelerometer bias' process noise, and construct your Kalman filter as usual.
Because you are describing the dependency of the observed object motion on the accelerometer bias in your state transition matrix, that dependency will make it into the covariance ($\mathbf P$) matrix. That will, in turn, affect the Kalman gain, which means that deviations in the object's predicted motion from measured will affect the observed accelerometer bias.