The problem is that my finite difference in the original post does not conserve mass. Mass conservation requires I use aconserving finite volume methoddifference schemes exist, rather than any arbitrary spatial discretizationbut I've found it easier to use and think in terms of finite volumes instead.
For exampledetails, see the finite volume textbook by Leveque
Below, I re-do the discretization, using finite volumes.
I can discretize the flux in the following way, considering the fluxes at the edges of each finite volume.
$$\frac{\partial \Phi}{\partial x}|_j \approx \frac{\Phi_{j+1/2} - \Phi_{j-1/2}}{\Delta x}$$
For a first-order upwind scheme, we can approximate the fluxes at the edges of the finite volumes in the following way
$$\Phi_{j+1/2} = a_j u_j$$ $$\Phi_{j-1/2} = a_{j-1} u_{j-1}$$
so,
$$\frac{\partial u_j}{\partial t} = - \frac{a_j u_j - a_{j-1} u_{j-1}}{\Delta x}$$
On the left boundary ($j=1$)
$$\frac{\partial u_1}{\partial t} = - \frac{a_1 u_1}{dx} + \frac{\Phi_{1/2}}{\Delta x}$$$$\frac{\partial u_1}{\partial t} = - \frac{a_1 u_1}{\Delta x} + \frac{\Phi_{1/2}}{\Delta x}$$
$\Phi_{1/2}$ is the flux at the left boundary. On the right boundary ($j = n$, where $n$ is the total number of finite volumes)
$$\frac{\partial u_n}{\partial t} = - \frac{\Phi_{n+1/2}}{dx} + \frac{a_{n-1} u_{n-1}}{\Delta x}$$$$\frac{\partial u_n}{\partial t} = - \frac{\Phi_{n+1/2}}{\Delta x} + \frac{a_{n-1} u_{n-1}}{\Delta x}$$
$\Phi_{n+1/2}$ is the flux at the right boundary.
This finite volume upwind scheme conserves mass to machine precision.