The dynamics of viscous immiscible pressure-driven multilayer flows in channels are investigated using a combination of modelling, analysis and numerical computations. More specifically the particular system of three stratified layers with two internal fluid-fluid interfaces is considered in detail in order to identify the nonlinear mechanisms involved due to multiple fluid surface interactions. The approach adopted is analytical/asymptotic and is valid for interfacial waves that are long compared to the channel height or individual undisturbed liquid layer thicknesses. This leads to a coupled system of fully nonlinear partial differential equations of Benney-type that contain a small slenderness parameter that cannot be scaled out of the problem. This system is in turn used to develop a consistent coupled system of weakly nonlinear evolution equations, and it is shown that this is possible only if the underlying base-flow and fluid parameters satisfy certain conditions that enable a synchronous Galilean transformation to be performed at leading order. Two distinct canonical cases (all terms in the equations are of the same order) are identified in the absence and presence of inertia, respectively. The resulting systems incorporate all the active physical mechanisms at Reynolds numbers that are not large, namely, nonlinearities, inertia-induced instabilities (at non-zero Reynolds number) and surface tension stabilisation of sufficiently short waves. The coupled system supports several instabilities that are not found in single long-wave equations including, transitional instabilities due to a change of type of the flux nonlinearity from hyperbolic to elliptic, kinematic instabilities due to the presence of complex eigenvalues in the linearised advection matrix leading to a resonance between the interfaces, and the possibility of long-wave instabilities induced by an interaction between the flux function of the system and the surface tension terms. All these instabilities are followed into the nonlinear regime by carrying out extensive numerical simulations using spectral methods on periodic domains. It is established that instabilities leading to coherent structures in the form of nonlinear travelling waves are possible even at zero Reynolds number, in contrast to single interface (two-layer) systems; in addition, even in parameter regimes where the flow is linearly stable, the coupling of the flux functions and their hyperbolic-elliptic transitions lead to coherent structures for initial disturbances above a threshold value. When inertia is present an additional short-wave instability enters and the systems become general coupled Kuramoto-Sivashinsky type equations. Extensive numerical experiments indicate a rich landscape of dynamical behaviour including nonlinear travelling waves, time-periodic travelling states and chaotic dynamics. It is also established that it is possible to regularise the chaotic dynamics into travelling wave pulses by enhancing the inertialess instabilities through the advective terms. ...