We explore the capability of lattice Boltzmann equation (LBE) method for complex fluid flows involving turbulence, mixing, and reaction. In the first study, LBE schemes for binary scalar mixing and multi-component reacting flow with reactions are developed. Simulations of initially non-premixed mixtures yield scalar probability distribution functions that are in good agreement with numerical data obtained from Navier-Stokes (NS) equation based computation. One-dimensional chemically-reacting flow simulation of a premixed mixture yields a flame speed that is consistent with experimentally determined value. The second study involves direct numerical simulation (DNS) and large-eddy simulation (LES) of decaying homogenous isotropic turbulence (HIT) with and without frame rotation. Three categories of simulations are performed: (i) LBE-DNS in both inertial and rotating frames; (ii) LBE-LES in inertial frame; (iii) Comparison of the LBE-LES vs. NS-LES. The LBE-DNS results of the decay exponents for kinetic energy k and dissipation rate ε, and the low wave-number scaling of the energy spectrum agree well with established classical results. The LBE-DNS also captures rotating turbulence physics. The LBE-LES accurately captures low-wave number scaling, energy decay and large scale structures. The comparisons indicate that the LBE-LES simulations preserve flow structures somewhat more accurately than the NS-LES counterpart. In the third study, we numerically investigate the near-field mixing features in low aspect-ratio (AR) rectangular turbulent jets (RTJ) using the LBE method. We use D3Q19 multiple-relaxation-time (MRT) LBE incorporating a subgrid Smagorinsky model for LES. Simulations of four jets which characterized by AR, exit velocity, and Reynolds number are performed. The investigated near-field behaviors include: (1) Decay of mean streamwise velocity (MSV) and inverse MSV; (2) Spanwise and lateral profiles of MSV; (3) Half-velocity width development and MSV contours; and (4) Streamwise turbulence intensity distribution and spanwise profiles of streamwise turbulence intensity. The computations are compared against experimental data and the agreement is good. We capture both unique features of RTJ: the saddle-back spanwise profile of MSV and axis-switching of long axis from spanwise to lateral direction. Overall, this work serves to establish the feasibility of the LBE method as a viable tool for computing mixing, combustion, and turbulence.