Este documento presenta un algoritmo para el modelado numérico de la dispersión de burbujas que ocurre en la capa cercana a la superficie del océano superior bajo la acción de ondas superficiales bidimensionales (2D) no rompientes. El algoritmo se basa en un enfoque Euleriano-Lagrangiano donde se resuelven las ecuaciones completas de Navier-Stokes en 3D para el flujo portador inducido por una superficie de agua ondulada en un marco Euleriano, y las trayectorias de las burbujas individuales se rastrean simultáneamente en un marco Lagrangiano, teniendo en cuenta el impacto de las burbujas en el flujo portador. Los diámetros de las burbujas se consideran en el rango de 200 a 400 micrones (por lo tanto, microburbujas), y se descuidan los efectos relacionados con la deformación y disolución de las burbujas en el agua. El algoritmo permite la evaluación de los perfiles instantáneos, así como estadísticamente estacionarios, promediados en fase del flujo turbulento del portador, la concentración de burbujas (fracción de vacío) y los flujos de fracción de vacío para diferentes regímenes de flujo, tanto con como sin deriva superficial inducida por el viento. Los resultados de las simulaciones muestran que las burbujas son capaces de aumentar la turbulencia del flujo portador, en comparación con el flujo sin burbujas, y que las fluctuaciones de la velocidad del agua vertical se incrementan principalmente, y cada vez más, por burbujas más grandes. Los resultados también muestran que la dinámica de las burbujas está gobernada por la flotabilidad, la fuerza de aceleración del fluido circundante y la fuerza de arrastre, mientras que el impacto de la fuerza de sustentación sigue siendo insignificante.