Introducción

La turbulencia sigue siendo uno de los problemas no resueltos más importantes de la física clásica (^{She, 1991}). Nuestro principal interés en la turbulencia se centra en la comprensión de la dispersión de contaminantes en atmósferas tropicales y su impacto en episodios de contaminación (Ren, et al., 2018). El transporte, dispersión y deposición de emisiones en el trópico están en gran parte determinados por convección húmeda y, frecuentemente, por fenómenos meteorológicos locales y regionales modulados por condiciones sinópticas que condicionan la cobertura nubosa y, por lo tanto, la radiación solar. La actividad humana en las ciudades se desarrolla mayoritariamente en la capa de superficie, en la cual ocurre el intercambio de energía, de momento, de gases y de aerosoles entre la superficie y la atmósfera. La actividad humana y el cambio del uso del suelo, que altera el albedo, la rugosidad aerodinámica, la altura del dosel y la homogeneidad de la superficie, genera flujos de calor suficientemente elevados como para modificar sustancialmente el balance energético (^{Roth, 2017}). Estas interacciones entre la superficie y la atmósfera producen movimientos turbulentos a distintas escalas espaciotemporales, que se manifiestan en remolinos intermitentes (She, 1991), coherentes (^{Starkenburg, et al., 2013}; ^{Christen, et al., 2007}), aniso-trópicos, o en ráfagas de viento.

Para estimar el transporte y dispersión de contaminantes en la atmósfera se emplean modelos eulerianos y lagrangianos, en los cuales los efectos de la turbulencia en la dispersión generalmente se han modelado con base en teorías de relación gradiente-transporte y de similitud (^{Arya, 1999}; ^{Hanna, 1979}; ^{Shao, 1992}; ^{Lin, et al., 2003}), entre las que se destaca la hipótesis de similitud de Monin-Obukhov (Arya, 1999). Las parametrizaciones de la turbulencia en estos modelos, y por tanto su habilidad para describir la dispersión de contaminantes, típicamente se ha basado en observaciones hechas en terrenos rurales, cuya superficie es mucho menos rugosa que la de las grandes ciudades (^{Christen, et al., 2007}), y se ha validado en latitudes medias y altas, pero muy poco en el trópico, y mucho menos en áreas urbanas tropicales de gran complejidad topográfica, sobre las cuales hay un enorme vacío de información.

Además, la comprensión de la turbulencia atmosférica también es importante en campos tan disímiles como la energía eólica (^{Hansen, et al., 2012}; ^{Barlas, et al., 2018}), la aviación y la ingeniería aeronáutica (^{Richardson, et al., 2014}), y la óptica atmosférica (^{Pinzón, et al., 2016}).

En la revisión bibliográfica realizada se pudo establecer que este es el primer estudio observacional directo, es decir, basado en mediciones de alta frecuencia, de la turbulencia en Colombia, y uno de los pocos llevados a cabo en ciudades tropicales y subtropicales (^{Yadav, et al., 1996}; ^{Roth, et al., 2017}) y, en general, en el trópico (^{Magnago, et al., 2010}; Silva da Santana, et al., 2015). En el único reporte previo sobre Colombia, ^{Montoya, et al. (2004)} estimaron varios parámetros micrometeorológicos para Bogotá (flujo de calor sensible, velocidad de fricción, temperatura de escala turbulenta y longitud de Obukhov), usando mediciones de baja frecuencia de la temperatura del aire a 2 m y 10 m de altura, y de la velocidad del viento a 10 m de altura en un periodo corto (3 días).

Las preguntas científicas que guían nuestra investigación son las siguientes: a) ¿cuáles son la magnitud, la variabilidad espaciotemporal y la distribución de frecuencia y espectral de la turbulencia atmosférica en Colombia? b) ¿qué factores meteorológicos, topográficos y de cobertura controlan esta variabilidad? c) ¿cuál es la estacionalidad y la variabilidad del intercambio entre superficie y atmosfera caracterizado como altura de capa de mezcla (^{Ortiz, 2016}) y flujos de momento, calor y masa, particularmente de gases de efecto invernadero, medidos mediante covarianza de remolinos? d) ¿qué tan apropiada es la parametrización de la turbulencia en los modelos de dispersión de contaminantes para las condiciones de Colombia y cómo se podría mejorar (Ortiz, 2016)? El presente análisis contribuye a responder la primera pregunta mediante la caracterización estadística y espectral de la turbulencia atmosférica en la capa de superficie en un sitio meteorológicamente representativo del área metropolitana de Bogotá.

Materiales y métodos

Sitio y período de medición. Bogotá es la ciudad más poblada de Colombia y una de las más pobladas de Suramérica, con una población estimada para el 2015 (año representativo para las mediciones reportadas) de 7.88 millones de habitantes (Mhab) en su área urbana y de 9.36 Mhab en su área metropolitana (Departamento Administrativo Nacional de Estadística - DANE, 2015). Se encuentra ubicada en la Sabana de Bogotá, una altiplanicie a ~2600 m de altitud, asentada en la trifurcación oriental de la cordillera de los Andes en Colombia (cordillera Oriental) (Figura 1SA-B, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3573). Como resultado de su posición en los Andes tropicales, su complejidad topográfica y la interacción entre la convección húmeda y efectos atmosféricos de escala sinóptica, mesoescala (vientos valle-montaña y anabáticos-catabáticos) y locales (canalización), la circulación atmosférica en Bogotá es relativamente compleja (^{Guerrero, 2013}; Guerrero & Jimenez, 2014; ^{Ortiz, et al., 2015}). La Sabana de Bogotá tiene un régimen bimodal de precipitación con picos climatológicos en abril y octubre y mínimos en enero y julio (^{Jaramillo-Robledo & Chaves-Córdoba, 2000}). Gran parte de la humedad que se precipita en la Sabana es transportada desde la cuenca del río Orinoco (Orinoquia), con muy pocos ciclos de evaporación-condensación según el análisis isotópico (^{Saylor, et al., 2009}). Las observaciones operacionales del Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM) indican que las nubes sobre Bogotá son mayoritariamente advectivas durante la primera temporada de lluvias (abril-mayo), es decir, formadas en la Orinoquia y transportadas hasta la Sabana, y en su mayoría convectivas durante la segunda temporada de lluvias (octubre-noviembre), es decir, generadas regionalmente (^{Arango, 2018}). El promedio de la precipitación acumulada anual en el Parque Metropolitano Simón Bolívar fue de 980 mm entre los años 2004 y 2013, siendo los meses de octubre y noviembre los de mayor precipitación mensual, con promedios mensuales superiores a los 120 mm acumulados. Las mayores velocidades del viento se presentan de julio a septiembre, cuando Bogotá está influenciada por los vientos alisios del este (Ortiz, et al., 2015).

Las mediciones de la turbulencia atmosférica se llevaron a cabo en la estación Parque Simón Bolivar (4°39'30.48" N, 74°5'2.38" O, 2577 m s.n.m.m.), conocida también como Centro de Alto Rendimiento (IDRD) de la Red de Monitoreo de la Calidad del Aire de Bogotá (RMCAB: http://ambientebogota.gov.co/red-de-calidad-del-aire). El Parque Simón Bolívar es una zona plana ubicada en el centro geográfico de Bogotá. Mediante un análisis de información mutua de las series temporales meteorológicas y de la calidad del aire de la RMCAB, y resultados de simulación meteorológica y lagrangiana estocástica, ^{Guerrero (2013)} y Guerrero & Jimenez (2012; 2014) encontraron que el Parque Simón Bolívar es uno de los sitios más representativos de la atmósfera de Bogotá, ya que está meteorológicamente influenciado por masas de aire que atraviesan la ciudad desde diferentes direcciones. Las mediciones aquí reportadas pueden considerarse altamente influenciadas por el dosel urbano, ya que el Parque está rodeado de edificaciones y su área verde (~4 km²) es muy pequeña comparada con la de Bogotá. El Parque puede dividirse en siete sectores según el uso del suelo (^{Morales-Rincon, et al., 2015a}) (Figura 1S-C, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3573). Los sectores 1 y 2 tienen coberturas homogéneas de baja rugosidad (canchas de futbol y campo de golf con pastos de menos de ~0.3 m). El sector 3 es ligeramente más heterogéneo y rugoso (pastos menos manejados y árboles dispersos). Los sectores 4, 5, 6 y 7 son los más heterogéneos y rugosos, ya que tienen pastos menos manejados y árboles dispersos entre canchas deportivas, edificios, algunos coliseos y una zona de parqueaderos. El fetch, es decir la distancia de cobertura homogénea viento arriba de la torre, para los sectores 1, 2, 3 y 6 es mayor a ~400 m, en tanto que para los sectores 4, 5 y 7 es mayor a 100 m (Figura 1S-C, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3573). La rugosidad del Parque se estima en ~0.5 m.

Las mediciones ("campaña" de medición) se llevaron a cabo entre el 25 de agosto y el 24 de noviembre del 2014. La ausencia de datos en el período comprendido entre el 24 de octubre y el 9 de noviembre de 2014 se debió a fallas en el suministro de potencia.

Instrumento de medición. El anemómetro sónico tridimensional (3D) de alta frecuencia (AS3D; Windmaster Pro 3D, Gill Instruments) instalado en el Parque hace parte de un sistema de medición de flujos de CO₂ y H₂O mediante covarianza de remolinos (eddy covariance, EC), el cual operó en Bogotá durante un período corto (Morales, et al., 2015a) antes de ser trasladado a un ecosistema agrícola (Morales, et al., 2015b; 2019). El sistema EC incluye el anemómetro AS3D, un analizador infrarrojo no dispersivo para la medición a alta frecuencia de la razón de mezcla de CO₂ y H₂O, y una estación meteorológica auxiliar (tem-peratura, humedad relativa de aire y suelo, precipitación, y radiación global y neta). El AS3D puede registrar datos hasta a 32 Hz, con una resolución nominal de 0.01 m s^-1, y una precisión (reportada como coeficiente de variación) de la velocidad escalar superior a una media cuadrática (root mean square, RMS) de 1.5 % en 10 s de medición a 12 m s^-1 (típicamente mejor que 1.3% RMS) Mc Kay, 2018, comunicación personal. A partir de esta información, se estima que la precisión normalizada (equivalente a ruido blanco) del AS3D es de 4.1 % Hz^-1/2, lo que equivale a una precisión de 13 % a 10 Hz, frecuencia de muestreo de los datos aquí reportados. Nuestra frecuencia máxima de observación según el teorema de Nyquist-Shannon (≤ 5 Hz) fue muy baja para poder observar la microescala de Kolmogorov (~10-30 Hz), pero suficientemente elevada para la microescala de Taylor (~0.01-5 Hz). El AS3D se instaló en la torre meteorológica del Parque a una altura no estándar de z = 7.8 m y al final de un brazo de 1 m de largo apuntado en la dirección NE, con el fin de garantizar que no afectara las mediciones del anemómetro propio de la estación (anemómetro de copa y veleta estándar, ACE; 034B Wind Sensor, Met One Instruments), instalado en la cima de la torre a 10 m de altura y sin obstáculos. El ACE tiene una exactitud nominal de ±4° y ±0.1 m s^-1 por encima del umbral de detección de 0.4 m s^-1. El AS3D registró mediciones a 10 Hz de cada uno de los componentes del viento, en tanto que el ACE registró la velocidad y la dirección horizontales del viento cada 10 minutos. Se compararon las rosas de los vientos de los dos anemómetros con una finalidad informativa, y se analizó la dependencia de las intensidades turbulentas de la dirección y la velocidad del viento medidas por el AS3D, con el fin de detectar posibles obstáculos.

Dado que los remolinos son transportados por el viento, la turbulencia observada por el AS3D (receptor) proviene, en general, de zonas aire arriba (upwind) de la torre meteorológica, cuyo mapa de influencia relativa en la medición se conoce como huella (footprint). En una investigación anterior (^{Morales-Rincon, et al., 2015a}), se calcularon y promediaron las huellas de observación del AS3D en cada uno de los cuatro meses de medición, las cuales se reproducen en el material suplementario (Figura 1S D-G). El cálculo se hizo en EddyPro (http://www.licor.com/ EddyPro) usando la parametrización de ^{Kljun, et al. (2004)}. En general, el "campo de visión" (field-of-view) del AS3D tuvo un alcance de ~600 m alrededor de la torre. El sector 3 fue el de mayor influencia durante agosto y septiembre, y los sectores 7, 1 y 2, los de mayor contribución durante octubre y noviembre.

Caracterización de la turbulencia. El AS3D midió los tres componentes ortogonales del viento, u, v, y w (este o zonal, norte o meridional, y hacia arriba en la vertical, respectivamente) cada 0.1 s. En terrenos suficientemente planos como el del Parque, se espera que para un tiempo de promediación apropiado, el promedio de la velocidad vertical del viento sea cero ( ). Sin embargo, es inevitable que haya pequeñas pendientes en el terreno y que errores menores de instalación impidan que el anemómetro sea instalado en una posición perfectamente vertical. Con el fin de evitar que esto genere sesgos en la velocidad vertical y flujos escalares artificiales, se aplicó rotación (tilt correction) de tres ejes (^{Foken, 2008}) cada 10 minutos (tiempo de promediación para la descomposición de Reynolds, t _avg ), previa eliminación de los picos de señal (despiking) (^{Starkenburg, et al., 2016}). Esta rotación estableció dos nuevos ejes ortogonales en la horizontal, uno en la dirección promedio del viento durante el período de promediación, llamado "longitudinal" u "horizontal", y el perpendicular a este, llamado "lateral", y un eje vertical corregido, ortogonal a los dos horizontales. La rotación preserva la notación de las variables: u, longitudinal, v, lateral, y w, vertical corregida (libre de sesgos, ). Los componentes advectivos ( ), se determinaron promediando las mediciones de alta frecuencia durante t _avg = 10 minutos. Los componentes turbulentos (u',v',w') se calcularon como los residuales alrededor de estos promedios usando las ecuaciones 1, 2 y 3 (por ejemplo, u' = u - ),

Inicialmente se hizo la descomposición de Reynolds con t _avg = 30 minutos, tiempo que resultó ser demasiado largo para la alta turbulencia lateral (^{Hanna, 1983}) en Bogotá, por lo que generaba artificialmente turbulencia elevada y picos espectrales de baja frecuencia, similares a los observados por ^{Magnago, et al. (2010)}. La magnitud y la variabilidad de los componentes turbulentos se expresaron como desviaciones estándar (sigmas) o medias cuadráticas, por ejemplo, en el caso de la velocidad longitudinal, Ecuación 4 donde u es la velocidad longitudinal promedio, y u _k y u’ _k son, respectivamente, la velocidad medida a alta frecuencia y su componente turbulento en el instante k-ésimo (a 10 Hz esto implica que para el cálculo de cada componente advectivo a t _avg = 10 min, se requieren n = 6000 mediciones). La magnitud del vector viento, (velocidad escalar), observada con el AS3D se calculó usando la Ecuación 5,

Las intensidades turbulentas por componente representan la relación entre la desviación estándar del componente turbulento y la velocidad escalar (comparable a coeficientes de variación), tal como se presenta en la ecuación 6 para el componente longitudinal,

La energía cinética turbulenta (turbulent kinetic energy, TKE) (ecuación 7) es la energía cinética media por unidad de masa disipada por los remolinos, cuya dependencia del tamaño de los remolinos se describe mediante distribuciones espectrales de energía (^{Arya, 1999}),

En el rango de frecuencia de interés para la micrometeorología (~30 minutos), la TKE de los remolinos se transfiere en cascada de las escalas grandes a las pequeñas. El ciclo de vida de los remolinos en la atmósfera incluye tres escalas de energía (^{Foken, 2008}; ^{Katul & Chu, 1998}): 1) producción (~10^-4 Hz), en la que las energías mecánica y convectiva, entre otras, forman y alimentan grandes remolinos en la atmósfera, por lo que son las propiedades de esta escala las que determinan la tasa de disipación de la energía (^{Mouri, et al., 2012}); 2) subrango de transferencia inercial de energía (microescala de Taylor, ~0,01-5 Hz), en la que los grandes remolinos se dividen progresivamente en pequeños remolinos, y 3) disipación viscosa (microescala de Kolmogorov, ~10-30 Hz), en la que los pequeños remolinos terminan por disiparse térmicamente vía transferencia de momento.

La función de autocorrelación, R, revela la similitud temporal de una variable mediante la correlación de su serie temporal con una copia rezagada de ella misma. La ecuación 8 presenta la autocorrelación para la turbulencia longitudinal en función del rezago de tiempo τ,

Por definición, R = 1 cuando τ = 0, valor que decae rápidamente a cero en el caso de la turbulencia. El tiempo de escala integral, T _i , es el tiempo máximo en el cual la variable aún se encuentra suficientemente correlacionada con sus valores anteriores, que, por definición, se calcula como el área bajo la curva de R (ecuación 9) (^{Arya, 1999}),

Cuando R(x) cambia a valores negativos, se usa la primera raíz de R como límite superior de integración (^{Swamy, et al., 1979}). Según la hipótesis de turbulencia estacionaria de Taylor (^{Arya, 1999}), el tiempo (T _i ) y la longitud de escala integral (L _i , ecuación 10) turbulentos representan la duración promedio de creación de los grandes remolinos y su longitud (~10-500 m) (^{Foken, 2008}),

El tiempo de la escala de Taylor, t _s (ecuación 11), es un estimador del tiempo con el cual los grandes remolinos se transforman en remolinos más pequeños en el subrango inercial,

y la longitud de la escala de Taylor, λ _s , representa el tamaño promedio de los remolinos más pequeños generados por la división de remolinos grandes (ecuación 12) (^{Arya, 1999}),

Por último, a partir de la transformada de Fourier, se determinaron los espectros de energía, S(n) (unidades [m²s^-1]), para los días de menor y mayor TKE diaria acumulada (equivalente a los días de menor y mayor TKE promedio), los cuales se redujeron a la densidad espectral adimensional (normalizada por la tasa de disipación de la energía y el inverso de la frecuencia), , y presentados en función de la frecuencia normalizada (adimensional), , donde n es la frecuencia [Hz], k es la constante adimensional de von Karman (k = 0.4), z es la altura de medición [m], y ε es la tasa promedio de disipación de la energía turbulenta por unidad de masa [m² s^-3] (ecuación 13),

donde L _e es la longitud de escala de la energía disipada o longitud característica de disipación [m], tradicionalmente definida como la longitud de correlación turbulenta de 2 puntos (^{Mouri, et al., 2012}), y C _ε es una contante adimensional de compleja definición (^{Therry & Lacarrere, 1983}), ya que depende de la configuración a gran escala del fluido. En condiciones de turbulencia estacionaria (^{Andreas, et al., 2008}), es geométricamente razonable estimar la longitud de escala de la energía disipada como el tamaño característico de los grandes remolinos, dejando implícita la constante adimensional, es decir como . Los espectros se promediaron para intervalos de 6 horas a partir de los componentes turbulentos derivados de la descomposición de Reynolds (10 minutos) y valores horarios de UJ y ε (Tabla 1S, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3577).

Discusión y análisis de resultados

Exposición del anemómetro sónico. La Figura 2S, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3574, muestra la rosa de los vientos observada por los dos anemómetros con un tiempo de promediación de 10 minutos. La mayor diferencia entre las frecuencias observadas por el AS3D y el ACE fue de 2 % en la dirección SSO. Las discrepancias para el resto de las direcciones no superaron el 1 % (0.8 % en promedio). La excelente comparabilidad a baja frecuencia del AS3D y el ACE, especialmente en la dirección SO (donde está ubicada la torre para el AS3D), muestra que la torre no constituye un obstáculo significativo para el AS3D. Más aun, un obstáculo significativo, incluida la torre misma, generaría estelas y turbulencia sustancial, lo que se manifestaría como una mayor intensidad turbulenta en la dirección del obstáculo. La gran isotropía y la homogeneidad de las intensidades turbulentas, que se presentan promediadas para toda la campaña en función de la dirección y la velocidad del viento horizontal (Figura 3S, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3575), es indicativa de la ausencia de obstáculos mayores en el sitio de medición. Estas dos evaluaciones confirman la buena exposición del AS3D. Las rosas muestran, además, que las direcciones predominantes del viento fueron E y SE. Las mediciones del ACE no se usaron en el resto del artículo.

Distribución de frecuencia y variabilidad temporal de los componentes turbulentos. Con el fin de evaluar que tan gaussiano es el comportamiento aleatorio de la turbulencia, las distribuciones de frecuencia empíricas de los componentes turbulentos instantáneos se normalizaron por la desviación estándar, por ejemplo, , y se compararon con la distribución normal (subíndice N: μ_Ν = 0, σ_Ν= 1) usando diagramas cuantil-cuantil. La Figura 4S, https://www.raccefyn.co/index.php/raccefyn/article/downloadSuppFile/697/3576, presenta los resultados para una selección aleatoria de 1000 intervalos de 10 minutos. Este análisis revela que la distribución gaussiana representa apropiadamente la turbulencia observada en las tres direcciones del viento en el rango ±1-sigma, pero con desviaciones significativas a partir de ±2-sigma. Estas desviaciones son, en promedio, simétricas tanto en signo como en magnitud. Sin embargo, no se espera que la turbulencia atmosférica sea puramente gaussiana, ya que los remolinos se correlacionan débilmente (^{She, 1991}) y de otro modo su autocorrelación sería cero para cualquier retraso de tiempo diferente a cero (^{Hanna, 1979}).

Las series temporales muestran la sucesión de ciclos diarios de turbulencia, expresada como desviaciones estándar y TKE de amplitud variable, pero típicamente caracterizados por máximos de 0.50.7 m s^-1 al mediodía local (12:00-14:00 tiempo local, LT), y mínimos de poco antes del amanecer (04:00-05:00 LT) (Figuras 1 y 2). Este análisis también revela menor turbulencia y variabilidad de la turbulencia en el periodo octubre a noviembre (meses de menor radiación) comparada con las registradas en septiembre (aunque fueron de buena de calidad, solo se hicieron mediciones válidas durante 6 días de agosto, por lo que no se incluyeron en la comparación). Resulta notable que la turbulencia lateral fuera comparable y, en algunos casos, ligeramente mayor que la turbulencia longitudinal. Las fluctuaciones de la dirección del viento en Bogotá a ~2 m fueron un factor ~4 mayores que las reportadas por ^{Hanna (1983)} para un terreno plano interrumpido por una pequeña colina en el sur de Idaho, Estados Unidos. Más aun, la similitud en forma funcional y la magnitud del perfil promedio diario (PPD) de los componentes turbulentos horizontales ( ) (Figura 2C-D), indican que la turbulencia fue esencialmente isotrópica en el plano horizontal. Por otro lado, los PPD de la turbulencia en la horizontal y en la vertical tienen formas similares, pero la magnitud de los componentes horizontales fue un factor ~2 mayor que en la vertical ( ). ^{Yadav, et al. (1996)} reportaron cocientes de en los intervalos 3-4 y 3-7 en condiciones estables e inestables, respectivamente, para un sitio de características comparables al Parque Simón Bolívar en una megaciudad subtropical (Nueva Delhi).

Figura 1 Series temporales de las observaciones meteorológicas realizadas en el Parque Simón Bolívar promediadas a 10 minutos. A) Precipitación; B) Radiación solar global; C) Temperatura del aire; D) Velocidad escalar del viento; E) TKE; F-H) Turbulencia expresada como desviaciones estándar de los tres componentes rotados de la velocidad del viento; I-K) Intensidad turbulenta por componente rotado. Las mediciones anemométricas (D-K) se hicieron con el AS3D, y las otras (A-C) con instrumentos auxiliares. 

Figura 2 Perfil promedio diario (PPD) de la turbulencia en el Parque Simón Bolívar durante los cuatro meses de medición (año 2014) expresada como B) TKE; C-E) desviaciones estándar para los tres componentes rotados de la velocidad del viento, y F-H) intensidad turbulenta por componente rotado. La sombra indica la variación 2-sigma alrededor del promedio. 

Magnitud del viento, intensidad turbulenta y energía cinética turbulenta. Durante la campaña, la velocidad escalar promediada a 10 minutos ( ) varió típicamente en el intervalo de 0.5-2.3 m s^-1 (0.01-5.06 m s^-1 para intervalo completo) (Figura 1D). Estas bajas velocidades del viento son comunes en el trópico (^{Yadav, et al., 1996}). El PPD de UJ muestra un máximo bien definido alrededor del mediodía y valores base entre la medianoche y el amanecer (Figura 2A). De forma similar, la TKE horaria varió típicamente en el intervalo de 0.1-2.0 m² s^-2 (0.002- 6.48 m² s^-2 para el intervalo completo) (Figuras 1E y 2B). Estas TKE son varias veces más altas que las observadas en la Amazonia brasileña bajo la influencia de corrientes de chorro nocturnas de bajo nivel (Silva da Santana, et al., 2015), y comparables con las reportadas para Nueva Delhi (Yadav, et al., 1996), donde las velocidades del viento a una altura menor de medición ( = 0.6-3.5 m s¹ a 4 m) fueron mayores que en Bogotá, lo cual indica que para los periodos reportados, la intensidad turbulenta (I) fue mayor en Bogotá que en Nueva Delhi. Por otro lado, el máximo diario de la TKE pasó de ser un valor meseta casi constante, entre ~10:00-15:00 LT, en septiembre a ser a un pico bien definido a las ~12:00 LT en noviembre (agosto se excluyó por las razones ya expuestas). Los valores mínimos nocturnos no presentaron variaciones significativas durante la campaña (Figura 2B).

I varió típicamente en los intervalos de 0.55 a 0.70 y de 0.20 a 0.45 para los componentes horizontales y vertical del viento, respectivamente (Figuras 1I-K y 2F-H). Estas intensidades turbulentas son aún más altas de lo esperado a las bajas velocidades del viento características de Bogotá. ^{Hansen, et al. (2012)} reportaron bajas I horizontales (~0.05-0.2 a ~5 m s^-1) a 20 m de altura en un parque eólico costa afuera en Dinamarca, lo cual es esperable sobre la superficie marina a velocidades altas. Para una colina aislada en terreno plano, ^{Hanna (1983)} reportó una dependencia hiperbólica de Iv con la velocidad del viento, con valores a ~2 m s^-1 en el intervalo de ~0.1-0.6 y una mediana de ~0.2, es decir, ~3 veces más baja que en Bogotá. Las intensidades turbulentas en los tres componentes reportadas más recientemente para dos sitios (urbano y parque suburbano) en Essen, Alemania (^{Weber & Kordowski, 2010}), son, en promedio, un factor ~2 más bajas que en Bogotá. Estas comparaciones sugerirían que la turbulencia "adicional" observada en Bogotá es atribuible a la turbulencia derivada de la convección húmeda, del dosel urbano y de la complejidad topográfica. La convección húmeda contribuiría a la turbulencia: 1) antes y durante la precipitación; 2) a través de la formación regional de nubes, particularmente en la segunda temporada de lluvias, y 3) de forma indirecta, debido al forzamiento rápido generado por la variación de la nubosidad y la radiación (^{Andreas, et al., 2008}). En la Figura 1 se presentan los períodos de precipitación (y su impacto en la radiación solar y la temperatura) durante la campaña, y sugiere un vínculo de estos con aumentos sustanciales de la turbulencia, particularmente la lateral. Sin embargo, identificar claramente la turbulencia debida a la convección húmeda es una tarea compleja (Andreas, et al., 2008) que excede los propósitos de este artículo.

El PPD para I _u y I _v fue casi constante a ~0.6, con valores ligeramente más elevados en el periodo de ~06:00 a 12:00 LT. El PPD para Iw fue ligeramente diferente, con un pico más marcado de ~0.45 a las ~09:00 LT, y valores más bajos (~0.2) antes de las ~06:00 que después de mediodía (~0.3). En los meses de mayor turbulencia (agosto y septiembre) durante la campaña, la turbulencia típicamente alcanzó su nivel máximo a media mañana (~10:00 LT), lo que podría estar asociado con la "ruptura" de la inversión térmica (^{Ortiz, 2016}), mientras que la velocidad del viento aumentaba hasta mediodía (~13:00 LT). Como consecuencia, el pico de intensidad turbulenta ocurre a media mañana (Figuras 2F-H).

Autocorrelaciones, tiempos y longitudes de escala. Las funciones de autocorrelación por componente (Ecuación 8) se calcularon para cada período de descomposición de Reynolds (10 minutos), y se promediaron por componente para cuatro períodos del día (00:00-06:00, 06:00-12:00, 12:00-18:00 y 18:00-24:00 LT), lo que permitió evidenciar que los componentes turbulentos tienden a durar correlacionados más tiempo en la noche que en el día (Figura 3). Los tiempos de la escala de Taylor variaron en los intervalos de 10 a 90 s y de 6 a 11 s para los componentes horizontales y vertical del viento, respectivamente. Los tiempos de escala integral para el componente vertical y los componentes horizontales variaron en los intervalos de 2 a 500 s y de 11 a 500 s, respectivamente (Figura 4A). Las longitudes de escala integral resultantes, asimilables al tamaño de remolinos en la zona de producción, variaron en los intervalos de 2 a 80 m y de 10 a 1000 m para el componente vertical y los componentes horizontales, respectivamente (Figura 4B).

Figura 3 Funciones de autocorrelación promedio para los componentes turbulentos: A) longitudinal, B) lateral, y C) vertical de la velocidad del viento en el Parque Simón Bolívar en los cuatro intervalos horarios del día 

Figura 4 Diagramas de cajas y bigotes para: A) los tiempos y B) las longitudes de escala integral y de Taylor para los tres componentes rotados de la turbulencia atmosférica en el Parque Simón Bolívar en los cuatro intervalos horarios del día durante la campaña. 

Por su parte, las longitudes de escala de Taylor, comparables con el tamaño de los remolinos más pequeños asociados a la transferencia inercial de energía, variaron en los intervalos de 1 a 30 m y de 5 a 80 m para el componente vertical y los componentes horizontales, respectivamente. Además, los tiempos de las escalas integral y de Taylor fueron significativa y sistemáticamente mayores en la noche que en el día. Cuando se considera la mayor velocidad del viento en el día, las longitudes de escala, asimilables al tamaño de los remolinos, resultaron ser mayores en el día que en la noche. Los remolinos más grandes ocurrieron en el período de 12:00 a 18:00 LT. En resumen, los grandes remolinos tuvieron duraciones (tamaños) modales de producción de ~11 s (~10 m) y ~160 s (~120 m) en la vertical y la horizontal, respectivamente. Por su parte, los remolinos pequeños, representativos del subrango inercial de energía, tuvieron tamaños (duraciones) modales de ~9 m (~9 s) y ~23 m (~24 s) en la vertical y la horizontal, respectivamente. Las longitudes de escala integral horizontales fueron de un orden de magnitud mayores que en la vertical, aproximadamente. Esta enorme disparidad de escalas es el resultado principalmente del confinamiento espacial de los remolinos, cuya relación de tamaños máximos hipotéticos sigue la proporción de longitud característica del área urbana a la altura de la capa de mezcla. En contraste, no hubo diferencias tan grandes entre los tiempos necesarios para transformar estos remolinos en remolinos más pequeños en la horizontal y la vertical y, por lo tanto, la diferencia entre las longitudes horizontales y vertical de la escala de Taylor fue mucho menor.

Tasa de disipación de energía turbulenta y espectros normalizados de energía. Nuestras estimaciones de la longitud de escala de la energía disipada son comparables a los valores reportados por ^{Therry & Lacarrère (1983)} para mediciones en la capa de superficie (L _ε /C _ε ~ 200 m para una altura de capa de mezcla de ~1000 m), y nuestras estimaciones resultantes de la tasa de disipación de energía por unidad de masa, de ~10^-3 m²/s³ y ~2^.10^-2 m²/s³, para los días de menor y mayor TKE, respectivamente (ver Tabla 1S), son comparables a los valores reportados por ^{Kaimal, et al. (1973)}.

La estabilidad atmosférica condiciona la forma funcional de los espectros normalizados al menos de tres formas (^{Kaimal, et al., 1976}): 1) cuanto mayor es la estabilidad menor es la turbulencia, lo que se refleja en un pico espectral menos pronunciado; 2) una mayor inestabilidad desplaza la frecuencia pico (frecuencia a la densidad espectral máxima) hacia frecuencias más bajas, y 3) los espectros en diferentes momentos del día y de los tres componentes turbulentos convergen, es decir se superponen, a frecuencias relativamente altas, independientemente de la estabilidad atmosférica, debido a que en el subrango inercial, donde ocurre la transferencia inercial de energía de los remolinos grandes a los pequeños, no hay forzamiento directo, es decir, la interacción de los remolinos con fuerzas mecánicas y convectivas es despreciable. La teoría de Kolmogorov predice que la densidad espectral normalizada en la zona inercial decrece dos décadas por cada tres décadas de aumento de la frecuencia normalizada, es decir, (Kaimal, et al., 1972; ^{Foken, 2008}).

La Figura 5 presenta los espectros normalizados de energía (suavizados mediante Kernel con una fracción de 0.75) para períodos de seis horas en los días de menor (10 de septiembre de 2014) y mayor TKE (30 de septiembre de 2014) durante la campaña. Los espectros de los cuatro periodos del día presentan la pendiente característica del subrango inercial, lo que marca su inicio, a frecuencias normalizadas de en la horizontal y de ׳ en la vertical. En el día de menor TKE, los espectros se superponen a partir de esas mismas frecuencias. Como resultado de una mayor variación de L _ε , los espectros son mucho menos compactos en el día de mayor TKE, en el cual tienden a converger solamente al final del subrango inercial y al comienzo de la microescala de Kolmogorov a .

Figura 5 Espectros reducidos de energía de los tres componentes rotados de la turbulencia atmosférica en el Parque Simón Bolívar en los cuatro intervalos horarios en los días de menor (09/10/2014) y mayor (30/09/2014) TKE durante la campaña 

Los picos espectrales caracterizan la zona de producción. Todos los espectros en la vertical tuvieron picos bien definidos a . En contraste, la mayoría de los espectros horizontales en el día de menor TKE son monotónicos. Para el día de mayor TKE, y los períodos alrededor de mediodía en el día de menor TKE, los picos ocurrieron a y . Las densidades espectrales y formas funcionales de los espectros son comparables a las ya clásicas reportadas por ^{Kaimal, et al. (1972)}, las cuales a su vez son muy similares en forma y magnitud a las reportadas para el norte de la India (^{Yadav, et al., 1996}) y el noreste del Brasil (^{Magnago, et al., 2010}). Para el día de mayor TKE, las frecuencias y densidades espectrales normalizadas para los componentes longitudinal y vertical corresponden a condiciones neutrales a ligeramente inestables en las curvas generalizadas de Kaimal, et al. (1972). En contrataste, para el componente lateral, incluso en el día de menor TKE, las frecuencias pico en Bogotá (_v’ = 0.01-0.07) son comparables a las del componente longitudinal, y un orden de magnitud menores que las reportadas en los tres artículos mencionados arriba (_v’ = 0.1-0.3), lo que concuerda con una turbulencia lateral elevada e isotropía horizontal. Además, las densidades espectrales son sustancialmente mayores, lo que caracterizaría los espectros laterales en Bogotá como significativamente inestables según Kaimal, et al. (1972).

Tal como era de esperar, las frecuencias pico se desplazaron hacia frecuencias menores en la zona de producción cuando la TKE aumentó, lo cual es indicativo de una mayor generación de turbulencia. Asimismo, independientemente de la TKE, las frecuencias pico aumentaron a medida que disminuía la producción de turbulencia entre el día y la noche. No se observaron diferencias sustanciales en la densidad espectral entre los días con menor y mayor TKE.

Conclusiones y perspectivas

Este es el primer estudio observacional directo y la primera caracterización estadística de la turbulencia en Colombia, y uno de los primeros sobre ciudades tropicales. La excelente comparabilidad a baja frecuencia con un anemómetro de copa estándar, y la isotropía polar de las intensidades turbulentas observadas a alta frecuencia, validan la buena exposición del anemómetro sónico y de las mediciones reportadas.

La distribución gaussiana de frecuencia representa apropiadamente la turbulencia observada en las tres direcciones del viento en el rango ±1-sigma.

Los ciclos diarios de la velocidad escalar, la turbulencia expresada como desviación estándar, y la TKE se caracterizan por máximos alrededor del mediodía y mínimos antes del amanecer. El presente estudio revela que, durante el periodo de medición, la turbulencia lateral en Bogotá fue elevada y comparable a la turbulencia longitudinal ( 0.2-1.3 m s^-1 PPD), lo que permite caracterizar su atmósfera como horizontalmente isotrópica. La turbulencia vertical ( 0.1-0.7 m s^-1 PPD) fue ~2 veces menor que la longitudinal, similar a lo reportado en la literatura. Las intensidades turbulentas observadas ( ) fueron ~2 a ~3 veces más altas que lo reportado para latitudes medias a las bajas velocidades del viento en Bogotá. Esta turbulencia "adicional" podría tener su origen en la convección húmeda, y en la complejidad topográfica y del dosel urbano.

Con base en análisis de autocorrelación, se estima que los remolinos grandes en la zona de producción, es decir al inicio de la cascada de energía, tuvieron longitudes características de ~120 m y ~10 m en la horizontal y la vertical, respectivamente.

El análisis espectral revela frecuencias y densidades espectrales de producción e inerciales concordantes con el análisis estadístico. El subrango inercial se extiende de frecuencias normalizadas de hasta de ( 100. La zona de producción está caracterizada por picos espectrales a . Las densidades espectrales y las formas funcionales de los espectros longitudinales y verticales son comparables a las reportadas en la literatura para condiciones neutrales a ligeramente inestables. En contraste, durante el periodo de medición, las frecuencias pico en Bogotá para el componente lateral (( _v’ = 0.01-0.07 ) fueron un orden de magnitud menores que las reportadas en la literatura, lo que concuerda con una turbulencia lateral elevada e isotropía horizontal.

Entre otros aportes, la caracterización hecha en este estudio provee las variables y parámetros necesarios para refinar y validar la simulación de huellas lagrangianas esto-cásticas en los Andes tropicales, con las cuales es posible, por ejemplo, ubicar fuentes de contaminación. Estas huellas se estiman a partir de trayectorias definidas cinemáticamente, por ejemplo, dx/ dt = u + u' para el componente zonal de una parcela de aire no sometida a deposición; uu generalmente se interpola a partir de campos meteorológicos simulados. En el caso más simple (turbulencia homogénea), el componente turbulento puede estimarse asumiendo "caminata aleatoria", , donde R _Lu es el coeficiente lagrangiano de autocorrelación (el cual puede estimarse a partir del tiempo de escala integral y la intensidad turbulenta), y u'' es la turbulencia puramente aleatoria (la cual puede estimarse a partir de la desviación estándar longitudinal y R _Lu ) (^{Lin, et al., 2003}; ^{Arya, 1999}).