Introducción

De acuerdo con ^{Grijalva-Contreras, Macías-Duarte, y Robles-Contreras (2011)}, los rendimientos anuales bajo invernadero de 24 híbridos de tomate en la región de Caborca, México, oscilaron entre 26.7 y 31.1 kg·m^-2, mientras que en el Reino Unido se reportan rendimientos de hasta 52.6 kg·m^-2 anuales (^{Antón et al., 2012}). Un elemento importante para elevar los rendimientos y la calidad de los productos agrícolas en invernadero es el control de las variables climatológicas, por lo que es necesario conocer el comportamiento de aquellos factores que las afectan directamente.

^{López-Cruz y Hernández-Larragoiti (2010)} han descrito modelos de caja negra que predicen de manera precisa variables como temperatura y humedad relativa dentro de un invernadero, en los cuales aplican modelos neurodifusos con R² entre 0.95 y 0.97. Asimismo, ^{Salazar-Moreno, López-Cruz, Rojano, Schmidt, y Dannehl (2015)} utilizaron un modelo dinámico de redes neuronales para predecir rendimientos de tomate en un invernadero semicerrado con buen desempeño, reflejado por una R² de 0.97. No obstante, este tipo de modelos tienen la desventaja de no incluir información sobre los procesos de transferencia de masa y energía, ni de los parámetros que intervienen en ellos; contrario a los modelos dinámicos del clima dentro del invernadero, que además proveen información sobre los procesos de intercambio de calor.

El balance de energía dentro del invernadero incluye todos los modos de transferencia de calor (por radiación térmica, conducción y convección), como se describe en los modelos empleados por ^{Al-Jamal (1994)}, ^{Arinze, Schoenau, y Besant (1984)}, ^{Baille (1999)}, ^{Boulard y Baille (1987)}, ^{Kindelan (1980)}, ^{Rodríguez, Berenguel, Guzman, y Ramírez-Arias (2015)}, ^{van Beveren, Bontsema, van Straten, y van Henten (2015a} y ^2015b) y ^{Walker (1965)}. El primer componente del balance de energía es la radiación solar que incide sobre la cubierta del invernadero, que puede ser transmitida, reflejada o absorbida. La proporción de radiación que atraviesa la cubierta se conoce como transmisividad, y depende de las características de la cubierta del invernadero y del tipo de radiación (directa o difusa) (^{Hernández, Escobar, & Castilla, 2001}). La cubierta aísla la atmósfera interna de las condiciones climáticas externas, por lo que funge como un enlace entre ambos entornos (^{Rodríguez et al., 2015}). Otro componente del balance energético es la ventilación (natural o forzada), la cual evita el calentamiento excesivo durante el día, asegura niveles mínimos de CO₂ y controla la humedad (^{Castilla-Prados, 2007}).

Por otro lado, la transpiración de las plantas produce una pérdida de calor en el invernadero, que depende de la concentración de vapor de agua, la conductancia de la transpiración, el índice de área foliar, la radiación neta del cultivo y la resistencia estomática que limita la transpiración (^{Bakker, Bot, Challa, & Braak, 1995}). La condensación del vapor de agua en invernadero, aunque no es muy grande, constituye otra pérdida de calor a considerar en el balance energético; así como el suelo, que constituye cerca de 10 % del total de pérdidas (^{Rodríguez et al., 2015}).

Al conocer las variables climatológicas externas (temperatura del aire, velocidad del viento y radiación solar), las propiedades de la cubierta, las especificaciones del cultivo y del invernadero, se puede plantear un conjunto de ecuaciones diferenciales y resolverlas de manera simultánea con el fin de describir el comportamiento de la temperatura del aire dentro del invernadero con respecto al tiempo.

Con base en lo anterior, los objetivos fueron: 1) desarrollar un modelo dinámico de balance de energía para predecir la temperatura del aire dentro de un invernadero experimental con cultivo de tomate, ubicado en la Universidad Autónoma Chapingo, México, b) obtener los valores óptimos de los principales parámetros involucrados en los procesos de transferencia de calor a través de la calibración del modelo y c) evaluar el modelo para su uso como una herramienta que permita optimizar el consumo de energía en sistemas de calefacción y controlar el clima en invernadero.

Materiales y métodos

Manejo del cultivo

El estudio se realizó en la Universidad Autónoma Chapingo, México, ubicada a los 19° 29’ latitud norte, 98° 53’ longitud oeste y 2,240 msnm, para lo cual se utilizó un invernadero de 8 m de ancho y 15 m de largo (120 m²), con cubierta de polietileno, malla anti insectos, dos ventanas laterales (6 m de ancho por 13 m de largo), una ventana cenital (15 m de largo por 1 m de ancho) y un sistema de apertura y cierre manual.

Se empleó el híbrido comercial de tomate (Solanum lycopersicum L.) “El Cid”, de hábito indeterminado y fruto tipo saladette. La siembra se realizó el 6 de marzo de 2016 en charolas de poliestireno de 200 cavidades con turba como sustrato, y el 24 de abril se trasplantó en sustrato de arena volcánica con densidad de 3.5 plantas·m^-2. Se fertilizó con la solución nutritiva Steiner (^{Steiner, 1961}) y el tutoreo se realizó de acuerdo con las necesidades del cultivo. La polinización se hizo en forma manual 15 días después del trasplante (ddt) y la cosecha inició 125 ddt.

Para la medición de las variables climatológicas se instalaron dos estaciones meteorológicas HOBO (Figura 1), una dentro del invernadero para registrar la temperatura interna, temperatura del suelo, radiación global y humedad relativa. Los sensores se ubicaron en el centro del invernadero y las mediciones se llevaron a cabo cada minuto. En el caso de la temperatura del techo, se instalaron dos sensores en la cubierta y los datos se almacenaron cada minuto en un sistema de adquisición de datos (datalogger). Además, se colocó un anemómetro en la parte exterior del invernadero a 2.5 m del suelo para medir la velocidad del viento.

Figura 1 Ubicación de las estaciones HOBO en el invernadero tipo cenital (8 m de ancho y 15 m de largo) con ventilación natural, localizado en la Universidad Autónoma Chapingo, México. 

Modelo dinámico de balance de energía

El modelo dinámico de balance de energía está basado en los modelos descritos por ^{de Zwart (1997)}, ^{van Beveren et al. (2015a} y ^2015b), ^{van Henten y Bontsema (2009)}, ^{van Ooteghem (2010)} y ^{Vanthoor, van Henten, Stanghellini, y de Visser (2011)}. En la Figura 2, se muestra la variable de estado temperatura del aire (Ta), las entradas (ganancia de calor debido a la radiación solar [Q _sol ] y ganancia de calor debido a la condensación de vapor de agua en el techo [Q _cond ]; ambas en W·m^-2) y salidas del sistema (pérdida de calor a través de la cubierta [Q _cub ], pérdida de calor debida a la ventilación [Q _vent ], pérdida de calor por la transpiración del cultivo [Q _trans ] y pérdida de calor a través del suelo [Q _suelo ]; todas en W·m^-2).

Figura 2 Flujos de energía que intervienen en un invernadero. 

Asimismo, en las Figuras 3 y 4 se presenta el comportamiento de las variables de entrada del modelo dinámico para 68 días (99,158 datos en total), correspondientes al periodo final de crecimiento del tomate (del 18 de junio al 26 de agosto de 2016).

Figura 3 a) Radiación solar (dentro), b) velocidad del viento (fuera) y c) humedad relativa (dentro), del 18 junio al 26 de agosto de 2016 (68 días). 

Figura 4 a) Temperatura dentro del invernadero, b) temperatura del suelo, c) temperatura de la cubierta y d) temperatura externa, del 18 junio al 26 de agosto de 2016 (68 días). 

El modelo de balance de energía en estado no estacionario para la temperatura del aire está dado por la siguiente ecuación diferencial ordinaria, donde la variable independiente es el tiempo (t) y la variable dependiente Ta:

VCpρdTadt=Qsol-Qcub-Qtrans-Qvent-Qsuelo+Qcond (1)

Si se consideran los flujos de energía por unidad de área (W·m^-2), la Ecuación (1) se modifica y en lugar del volumen (V, m³) se tiene h:

dTadt=1hCpρQsol-Qcub-Qtrans-Qvent-Qsuelo+Qcond (2)

donde h=VAs. La descripción de los parámetros y las variables incluidas en el modelo se muestran en los Cuadros 1 y 2.

Cuadro 1 Valor de los parámetros utilizados en la simulación. 

Símbolo	Descripción	Unidades	Valor
V	Volumen del invernadero	m3	630.00
As	Área del piso del invernadero	m2	120.00
Ac	Área de la cubierta	m2	345.3
Avent	Área de ventilación	m2	171.00
l0	Longitud de las ventanas	m	3.20
w0	Ancho de las ventanas	m	0.80
ps	Profundidad a la que se estima la temperatura del suelo	m	0.05
Cp	Calor especifico del aire	J·kg-1·oC-1	1000.00
ρ	Densidad del aire	kg·m-3	1.22
αc	Coeficiente de transferencia de calor de la cubierta	W·m-2·°C-1	5
τtot	Coeficiente de transmisión de la radiación solar	-	1
ks	Coeficiente de intercambio de calor a través del suelo	W·m-1·°C-1	5.75
k	Parámetro de transpiración específica del cultivo	-	0.40
kpo	Coeficiente de transferencia de masa	m·s-1	0.002532
rb	Resistencia aerodinámica	s·m-1	200
λ	Calor latente de vaporización	KJ·kg-1	2,500
IAF	Índice de área foliar	cm2·cm-2	2.2
P	Presión atmosférica	m·s-2	98.1
F 0	Coeficiente de pérdida de presión	-
Cf	Coeficiente de infiltración	m3·s-1·m-2	0.008361
Cd	Coeficiente de descarga	-	0.6497
Cw	Coeficiente de efecto del viento	-	0.14

Cuadro 2 Variables del balance de energía. 

Símbolo	Descripción	Unidades
Irad	Radiación solar externa	W·m-2
Qsol	Ganancia de calor debido a la radiación solar	W·m-2
Qcond	Ganancia de calor debido a la condensación	W·m-2
Qcub	Pérdida de calor a través de la cubierta	W·m-2
Qtrans	Pérdida de calor por la transpiración del cultivo	W·m-2
Qsuelo	Pérdida de calor a través del suelo	W·m-2
Qvent	Pérdida de calor debido a la ventilación	W·m-2
φcond	Conductancia de condensación	kg·s-1
Rn	Radiación neta a nivel cultivo	W·m-2
Ta	Temperatura interior	°C
Tc	Temperatura de la cubierta	°C
Ts	Temperatura del suelo	°C
Te	Temperatura exterior	°C
HRe	Humedad relativa externa	%
HRin	Humedad relativa interna	%
Xcul	Concentración de vapor de agua absoluto en el cultivo	kg·m-3
X	Concentración de vapor de agua en el invernadero	kg·m-3
Xsat	Concentración de vapor de agua saturado	kg·m-3
ɛ	Relación entre calor latente sensible	-
GE	Conductancia de la transpiración	m·s-1
ϕcond	Flujo de condensación	kg·s-1
rs	Resistencia estomática	s·m-1
ve	Velocidad del aire exterior	m·s -1
Tsv	Tasa de infiltración por unidad de área de piso	m3·m-2·s-1
Tv	Tasa de ventilación por unidad de área de piso	m3·m-2·s-1
t	Tiempo	s

Por otro lado, Q _sol se puede expresar como:

Qsol=τtotIrad (3)

Dado que la radiación solar se midió dentro del invernadero, se tomó el valor de 1 para el coeficiente de transmisión de la radiación solar (τ _tot ). Mientras que, la Q _cub está basada en el modelo de ^{van Beveren et al. (2015b)} (Ecuación 4).

Qcub=αcAcAsTa-Te (4)

Para el cálculo de la Q _trans de tomate en invernadero se utilizó el modelo propuesto por ^{Stanghellini y Jong (1995)}:

Qtrans=GEλχcul-χ (5)

donde se consideró el valor del calor latente de vaporización (λ) reportado por ^{Monteith y Unsworth (2013)} (λ = 2,500 J·g^-1). Mientras que, la conductancia de transpiración (G _E ) (^{van Beveren et al., 2015b}), la relación entre calor latente y sensible del contenido del aire saturado (ε), el parámetro de resistencia estomática (r _s ) y la radiación neta (R _n ) a nivel del cultivo (^{Bontsema et al., 2007}) están dadas por las siguientes ecuaciones:

GE=2IAF1+εrb+rs (6)

ԑ=0.7584e0.0518Ta (7)

rs=82+570e-kRnIAF1+0.023Ta-202 (8)

Rn=0.861-exp-0.7IAFQsol (9)

La concentración de vapor de agua en el cultivo (X _cul ) y la concentración del vapor saturado (X _sat ) se muestran en las Ecuaciones 10 y 11, respectivamente:

χcul=χsat+εrb2IAFRnλ (10)

χsat=5.563e0.0572Ta (11)

Por otro lado, la Q _vent se obtuvo con la ecuación descrita por ^{Ruiz-García, López-Cruz, Arteaga-Ramírez, y Ramírez-Arias (2015)}:

Qvent=Tv*ρ*Cp*Ta-Te    (12)

en donde la tasa de ventilación (T _v ; ^{Baeza et al., 2014}) y el coeficiente de descarga (C _d ; ^{Ruiz-García et al., 2015}) están dados por:

Tv=Avent2CdCw*ve+Tsv (13)

Cd=F0-0.5; F0=1.75+0.7e-[l0w0]/32.5 (14)

Sustituyendo los valores de ancho y largo de las ventanas, se obtiene un valor de C _d = 0.6497. En cuanto al coeficiente de efecto del viento (C _w ), ^{Ruiz-García et al. (2015)} reportan un valor de 0.16 para un invernadero con una superficie de 2,565 m²; mientras que ^{Roy, Boulard, Kittas, y Wang (2002)} y ^{Valera, Molina, y Álvarez (2008)} presentan un valor de 0.14 para un invernadero de 179 m², el cual se utilizó en este trabajo.

La primera parte de la Ecuación 13 es la ventilación por las ventanas, mientras que la segunda parte es la ventilación debida a las fugas de la construcción del invernadero, o tasa de infiltración (Ts _v ):

Tsv=Cf*As (15)

La ganancia de calor debida a la condensación se obtuvo con la ecuación de ^{Speetjens, Stigter, y van Straten (2010)}:

Qcond=λφcond    (16)

donde la conductancia de condensación (φ _cond ) se calculó de la siguiente manera:

φcond=kpo*ρ*Ac*χ-χ sat    (17)

Para esta ecuación, el coeficiente de transferencia de masa (k _po ) determina la tasa de condensación, y la X _sat y la concentración de vapor de agua en el invernadero (X) están dadas por:

χsat(Tc)=0.611[1+1.414sin⁡5.82e-3Tc]8.827 (18)

χ=0.611HRin χsat(Ta)100P (19)

Para la conversión de la humedad relativa a absoluta se utilizó la expresión propuesta por ^{Blasco, Martínez, Herrero, y Ramos (2007)} (Ecuación 20), la cual es una función que asegura que la humedad relativa nunca será mayor que 100 %.

HRin(Ta)= 100,  HRin>100HRin,  HRin≤100                    (20)

Finalmente, el Q _suelo se estimó con la ecuación de ^{van Ooteghem (2010)}. La profundidad a la que se estimó la temperatura del suelo fue a 0.05 m en la primera capa del suelo:

Qsuelo=ksAsTa-Tsps (21)

Para generar el modelo computacional se utilizó el ambiente de programación MATLAB/SIMULINK. Este consta de un programa principal donde se definen los parámetros (variables de entrada y salidas, así como las opciones de simulación) y un archivo C-MEX (S-function) donde se implementan las ecuaciones del modelo, que es llamado por un archivo mdl de SIMULINK. La simulación numérica se realizó con el método de integración tipo Runge-Kutta de 4º orden, conocido como método Dormand-Prince función ode45 de Matlab, con tamaño de paso de integración variable, y tolerancia relativa de 1x10^-8 y absoluta de 1x10^-10. Se dispuso de 99,158 datos en total para construir, calibrar y evaluar el modelo.

Resultados

Las variables climáticas, tanto internas (temperaturas del aire, suelo y techo, radiación solar, humedad relativa) como externas (temperatura y velocidad del viento), se obtuvieron cada minuto. De un total de 99,158 datos disponibles, 48 % (47,997 datos, del 18 junio al 22 julio de 2016) se utilizaron para construir y calibrar el modelo, y 52 % (51,153 datos, del 22 julio al 26 agosto de 2016) para evaluarlo.

Los resultados de la simulación del modelo para la temperatura del aire dentro del invernadero se compararon con los datos observados (Figura 5). Además, se obtuvieron diferentes medidas estadísticas de ajuste (Cuadro 3) para determinar la cercanía entre las mediciones de la temperatura del aire y las predicciones del modelo dinámico.

Figura 5 Comparación entre la temperatura dentro del invernadero medida y la simulada por el modelo antes de la calibración. 

En el modelo desarrollado, ocho parámetros están relacionados con las dimensiones del invernadero (V, A _s , A _c , A _vent , l ₀ , w ₀ , p _s , F ₀ ) y cuatro son constantes (C _p , ρ, λ, P); además, el τ _tot fue igual a 1, ya que la radiación solar se midió dentro del invernadero. Debido a lo anterior, no se realizó el análisis de sensibilidad, por lo que se hizo directamente la calibración del modelo para mejorar su desempeño con los parámetros restantes: resistencia aerodinámica (r _b ), coeficiente de transferencia de calor de la cubierta (α _c ), coeficiente de intercambio de calor a través del suelo (k _s ), parámetro de transpiración específica del cultivo (k), coeficiente de infiltración (C _f ) y coeficiente de efecto del viento (C _w ). Los valores iniciales de estos parámetros se obtuvieron de la literatura (Cuadro 4), los cuales se hicieron variar en diferentes porcentajes respecto a los valores nominales. El rango de variación que generó mejores resultados se muestra en el Cuadro 4.

El desempeño del modelo mejoró sustancialmente después de la calibración (Figura 6 y Cuadro 4), donde el porcentaje de cambio de la disminución del ECM, EAM, RECM y REAM fue de 70.32, 46.89, 45.47 y 13.23 %, respectivamente, y el cambio porcentual de incremento en la eficiencia fue 35.44 %.

Figura 6 Comparación entre la temperatura del aire dentro del invernadero medida y la simulada por el modelo después de la calibración. 

Una vez calibrado, el modelo desarrollado se puede utilizar como una caja negra para probar su desempeño con diferente información; es decir, también predice las respuestas cuando se alimenta de información nueva. A esto se le conoce como evaluación del modelo, la cual se implementó con el 52 % de los datos (Cuadro 4 y Figura 7). Por otro lado, el modelo de balance de energía predice en forma aceptable la variación de la temperatura del aire dentro del invernadero.

Figura 7 Predicción de la temperatura del aire con los valores de los parámetros estimados en la calibración y un nuevo conjunto de datos (evaluación del modelo). 

Discusión

En relación con los resultados de la calibración del modelo, el valor de los parámetros (Cuadro 3) se comparó con resultados obtenidos por otros investigadores. El valor de la r _b encontrado es mayor que el típico sugerido por ^{Bontsema et al. (2007)} para un invernadero holandés con cultivo de tomate (r _b = 200 s·m^-1), pero es similar al obtenido por ^{Qiu et al. (2010)} para tomate en un invernadero solar de China (r _b = 295 s·m^-1).

El C _w depende de las dimensiones del invernadero y de la velocidad del viento, por lo que después de la calibración se obtuvo un valor ligeramente superior (C _w = 0.2) al encontrado por ^{Ruiz-García et al. (2015)} (C _w = 0.16) para un invernadero de polietileno con cultivo de tomate. Sin embargo, se encuentra dentro del rango de valores reportados por ^{Molina-Aiz, Valera, Peña, Gil, y López (2009)} (entre 0.16 y 0.82). Asimismo, el C _f fue mayor (0.0125 m³·s^-1·m^-2) al obtenido por ^{Ruiz-García et al. (2015)}, lo que significa que el modelo sobreestima la tasa de infiltración.

Por otro lado, el α _c es superior (10 W·m^-2·°C^-1) al reportado por la Organización de las Naciones Unidas para la Alimentación y la Agricultura (^{FAO, 2002}) para invernaderos de plástico (entre 6 y 8 W·m^-2·°C^-1) y de cristal (8.8 W·m^-2·°C^-1 ). Mientras que el valor óptimo encontrado para el k _s (1 W·m^-1·°C^-1) es igual al reportado por ^{Valera et al. (2008)} para suelo seco, y se encuentra dentro del intervalo de los valores obtenidos por ^{Abu-Hamdeh y Reeder (2000)} para suelo arenoso (entre 0.58 y 1.94 W·m^-1·°C^-1) y suelo franco arenoso (entre 0.19 y 1.12 W·m^-1·°C^-1). Finalmente, el k fue menor al propuesto por ^{van Beveren et al. (2015b)}, lo que implica que se requiere mayor investigación para determinar el papel que cumple este parámetro, al usar un valor cercano a los nominales (reportados en la literatura), sobre la calidad de la predicción del modelo.

La sub y sobreestimación de los parámetros se puede atribuir a que en la construcción, calibración y evaluación del modelo, el área de ventilación se consideró constante, ya que el invernadero estuvo abierto la mayor parte del tiempo, únicamente se cerraban las ventilas en tiempo de lluvia de las 7:00 pm a 8:00-8:30 am. Esto último no se consideró en el modelo, y su inclusión podría mejorar su desempeño. Además, el índice de área foliar se consideró constante, pero su inclusión como variable de entrada podría mejorar la estimación de la transpiración.

En cuanto al planteamiento del modelo como una sola ecuación de balance de energía para estimar la temperatura interna en el invernadero, se considera que se puede mejorar al incluir balances en estado no estacionario a nivel de la cubierta, en la planta y en el suelo, y utilizar un sistema de ecuaciones diferenciales ordinarias y no solo una ecuación en el balance. Es decir, llevar a cabo balances en varias capas dentro del ambiente del invernadero (^{de Zwart, 1997}). Asimismo, se puede mejorar la calidad de predicción de la temperatura del aire dentro del invernadero si se considera el balance de materia asociado con la humedad del aire. Sin embargo, si se genera un modelo simple (con pocas variables de estado), este se puede usar fácilmente en la optimización y control del ambiente del invernadero (^{Seginer, van Straten, & van Beveren, 2017a} y ^2017b; ^{van Beveren et al., 2015a} y ^2015b).

Conclusiones

Se desarrolló un modelo de balance de energía para predecir la temperatura del aire dentro del invernadero y se identificaron los principales parámetros que influyen en el comportamiento térmico del invernadero, los cuales se calibraron. Los resultados obtenidos después de la calibración mejoraron de manera sustancial el desempeño del modelo. Las medidas de ajuste indican un modelo con eficiencias de 89 y 86 % después de la calibración y en la evaluación del modelo, respectivamente.

A pesar de que el parámetro de transpiración específica del cultivo (k) no está muy cercano a los valores reportados en la literatura, el modelo, en general, tiene un desempeño adecuado y puede utilizarse como una herramienta de apoyo para analizar el comportamiento de la temperatura dentro de un invernadero bajo diferentes condiciones en la región central de México.