Euler-Euler large eddy simulation (EELES) of transient dispersed turbulent bubbly flows in a laboratory scale square cross-sectioned bubble column has been presented. The main objective is to investigate the scale-adaptive of EELES model for predicting the dispersed bubbly flows. The sub-grid modeling is based on the Smagorinsky kernel with dynamic CS constant. The bubble induced turbulence and various interfacial forces including drag, lift, and virtual mass forces are incorporated in the model. Good quantitative agreement with previous experimental measurements by a two-camera particle image velocimetry (PIV) is obtained both for the fluctuating velocities and the mean velocities. More instantaneous details of two-phase flow characteristic in the bubble column have been captured by LES, including the multi-scale vortex structures, fluctuation characteristics of liquid motion and rotating buoyancy motion of bubbles. Predicted turbulence effective viscosity by unsteady Reynolds-Averaged Navier-Stokes (RANS) approach is nearly thirty times larger than LES, the turbulent dissipation is so strong that it cannot be used to predict the multi-scale fluctuation characteristics inside the bubble column. The classical -5/3 law of power spectrum densities (PSD) of liquid vertical velocity is hold properly in the low frequency region. For high frequency region, the classical -5/3 law does not hold properly, the decay is too fast, where the slope is around -25/3. The current EELES can be effectively used for studying the transient two-phase flow when a grid comparable to the bubble diameter is used, not just for the condition that the grid size is larger than the bubble diameter.