A series of experiments on turbulent bubbly channel flows observed bubble clusters near the wall which can change large-scale flow structures. To gain insights into clustering mechanisms, we study the interaction of a pair of spherical bubbles rising in a vertical channel through combined experiments and modelling. Experimental imaging identifies that pairwise bubbles of 1.0 mm diameter take two preferential configurations depending on their mutual distance: side-by-side positions for a short distance ($S<5$) and nearly in-line, oblique positions for a long distance ($S>5$), where $S$ is the mutual distance normalized by the bubble radius. In the model, we formulate the motions of pairwise bubbles rising at $Re=O(100)$, where $Re$ is the Reynolds number defined based on the bubble diameter. Analytical drag and lift, and semi-empirical spatio-temporal stochastic forcing, are employed to represent the mean acceleration and the fluctuation due to turbulent agitation, respectively. The model is validated against the experiment through comparing Lagrangian statistics of the bubbles. Simulations using this model identify two distinct time scales of the interaction dynamics which elucidate the preferential configurations. For pairs initially in line, the trailing bubble rapidly escapes from the viscous wake of the leading bubble to take the oblique position. Outside of the wake, the trailing bubble travels on a curved path with a slower velocity driven by potential interaction and horizontally approaches the leading bubble to become side by side. Moreover, statistical analysis identifies that the combination of the wake and the agitation can significantly accelerate the side-by-side clustering of in-line pairs. These results indicate positive contributions of liquid viscosity and turbulence to the formation of bubble clusters.