The transient motion of bubbly flows in vertical channels is studied, using direct numerical simulation (DNS) in which every continuum length and time scale is resolved. A simulation of a large number of bubbles of different sizes at a friction Reynolds number of 500 shows that small bubbles quickly migrate to the wall, but the bulk flow takes much longer to adjust to the new bubble distribution. Simulations of much smaller laminar systems with several spherical bubbles have been used to examine the full transient motion; those show a nonmonotonic evolution where all the bubbles first move toward the walls, and the liquid then slowly slows down, eventually allowing some bubbles to return to the center of the channel. Unlike the statistically steady state, where the flow structure is relatively simple and in some cases depends only on the sign of the bubble lift coefficient, the transient evolution is more sensitive to the governing parameters. Early efforts to use DNS results to provide values for the unresolved closure terms in a simple average model for the flow found by statistical learning from the data using neural networks are discussed. The prospect for using the results from simulations of large systems with bubbles of different sizes in turbulent flows for large eddy–like simulations are explored, including the simplification of the interface structure by filtering. Finally, preliminary results for flows undergoing topology changes are shown.