\section{Introduction}
\label{sec:intro}

A large amount of training data is often essential for neural model performance, especially, for large networks. Having more training data can help reduce overfitting and improve model robustness. However, preparing a large amount of annotated data is usually costly, labor intensive and time-consuming. Data augmentation \cite{Simard1998} is a useful technique for synthetic data generation, which is widely used in computer vision \citep{fawzi2016adaptive,d2017bridging,wang2019survey} and speech \citep{schluter2015exploring,ko2017study}.

However, due to the complexity of language, it is more challenging to apply data augmentation techniques to natural language processing (NLP). Unlike computer vision and speech, where handcrafted rules (such as rotation, cropping, masking, etc.) can be easily applied to transform original data, it is difficult to generalize such rules for languages. Although simple distortion usually does not change the semantics of visual information, deleting or replacing a single word could completely change the meaning of the sentence.

\begin{figure}[t!]
\centering
\includegraphics[scale=0.34]{image/linearization.png}
\caption{An example of labeled sentence linearization. All words and their tags are paired up by inserting tags before (or after) the words (\emph{O} tags removed).}
\label{fig:linearization}
\end{figure}

One successful method for data augmentation in NLP is \emph{back translation} \citep{sennrich-etal-2016-improving,fadaee2017data,dong-etal-2017-learning,wei2018fast}, where a translation model is used to translate monolingual sentences from target language to source language to generate synthetic parallel sentences. Other successful methods include: systematically reordering the dependents of some nodes in gold data to generate synthetic data for dependency parsing \citep{wang2016galactic}, leveraging knowledge base for question generation \citep{serban2016generating} and using simulation-based approach to generate a set of prerequisite toy tasks for QA \citep{weston2015towards}. Besides, synonym replacement, random deletion/swap/insertion, generation with VAE or pre-trained language models are also used in some NLP tasks \citep{kobayashi2018contextual,wei2019eda,anaby2020not,raille2020fast,kumar2020data}, but mainly for translation and classification tasks.

Compared with the above-mentioned downstream tasks like translation and classification, sequence tagging is more fragile when it is confronted with data augmentation noises due to the finer granularity of the (token-level) task. Annotating unlabeled data with a weak tagger, leveraging aligned bilingual corpora to induce annotation and synonym replacement are three attempted data augmentation methods for sequence tagging \citep{shang2018automated,yarowsky2001inducing,mathew2019biomedical}. Weakly labeled data will inevitably introduce more noise. Note that annotating unlabeled data with a weak tagger requires in-domain data and in-domain knowledge, otherwise it may suffer from domain-shift problem \cite{bari2020multimix}. Leveraging aligned bilingual corpora requires additional resources, which may not be available for low resource languages. Synonym replacement often relies on additional knowledge, e.g., WordNet \citep{miller1995wordnet}, which is a manually designed dictionary that may have low coverage (or not available) for low-resource languages.

In this work, we investigate data augmentation with a generation approach for sequence tagging tasks. We first linearize the labeled sentences as shown with an example in Figure~\ref{fig:linearization}. Then a language model (LM) is trained on the linearized data and used to generate synthetic labeled data. Unlike employing weak taggers to label unseen data, our method unifies the processes of sentence generation and labeling using a LM.
Concretely, a word and its tag in a pair (e.g., ``B-PER Jose'') are trained to be generated together, in which the tag-word pairs with high probability will be chosen by the LM during generation (\cref{subsec:lm_and_gen}). Our method does not require additional resources like WordNet. Nevertheless, if unlabeled data or knowledge bases are available, our method is also flexible to utilize these resources with a simple but effective conditional generation technique (\cref{subsec:conditional_generation}).

Although some recent work \citep{anaby2020not,raille2020fast,kumar2020data} also leverages LM for data augmentation, their methods are conditioned on sentence-level tags to generate or modify training data, hence applicable to classification tasks exclusively. To the best of our knowledge, we are the first to utilize generative language models to generate fine-grained synthetic data from scratch for sequence tagging tasks, which introduces a new paradigm for data augmentation in NLP. Furthermore, our method does not rely on large pre-trained models and in fact, it employs a simple one-layer recurrent language model \citep{sundermeyer2012lstm}, which is more convenient to train. Our method demonstrates  encouraging performance when trained on just a few thousand sentences (\cref{sec:exp}).

To verify the effectiveness of our method, we conduct extensive experiments on different sequence tagging tasks, including named entity recognition (NER), part-of-speech (POS) and end-to-end target based sentiment analysis (E2E-TBSA). Our method consistently outperforms the baseline methods in both supervised and semi-supervised settings. Different from the baseline methods, our method generates novel synthetic data from scratch, and thus introduces more diversity to reduce overfitting. For the semi-supervised settings, our method demonstrates strong ability to exploit useful information from unlabeled data and knowledge base.

\section{Background}
\label{sec:background}

\paragraph{Named Entity Recognition (NER)}
Named entities refer to phrases that are names of persons, organizations and locations, etc. in text. For example, ``\emph{[ORG U.N.] official [PER Ekeus] heads for [LOC Baghdad]} ".
Named entity recognition is an important task of information extraction and it aims to locate and classify named entities in text into the predefined types \citep{mikheev1999named,sang2003introduction,li2020survey}. It is a challenging task for two reasons  \citep{lample2016neural}: 1) in most languages and domains, the amount of manually labeled training data for NER is limited; 2) it is difficult to generalize from this small sample of training data due to the constraints on the kinds of words that can be names.

\paragraph{Part-of-Speech (POS) Tagging}
Part-of-speech tagging consists of assigning a tag that represents a grammatical class to each word in a given sentence. It is a critical component of most NLP systems and is fundamental to facilitate downstream tasks such as syntactic parsing \citep{schutze1993part} and opinion analysis \citep{liu-etal-2015-fine}.
The current state-of-the-art POS taggers can achieve over 97.80\% accuracy on PTB-WSJ \citep{akbik2018contextual,bohnet2018morphosyntactic} and yield over 96.50\% average test accuracy across 21 high-resource languages in UD 1.2 \citep{heinzerling2019sequence}. However, one of the problems with current POS taggers is that their accuracy can decrease significantly on low-resource languages and rare words \citep{plank-etal-2016-multilingual,yasunaga2018robust}.

\paragraph{Target Based Sentiment Analysis}
The target based sentiment analysis is a fundamental task of sentiment analysis and it aims to detect the opinion targets in sentences and predict the sentiment polarities over the targets \citep{liu-etal-2015-fine,chen-EtAl:2017:EMNLP20171,li2018aspect,li2019unified}.  For example, \emph{``\textbf{USB3 Peripherals} are noticeably less expensive than the \textbf{ThunderBolt ones}"}. In this sentence, two opinion targets were mentioned, namely ``\emph{\textbf{USB3 Peripherals}}" and ``\emph{\textbf{ThunderBolt ones}}" and the user expresses a positive sentiment over the first, and a negative sentiment over the second. \citet{li2019unified,li2019bert-e2e-absa} propose an end-to-end solution (E2E-TBSA) of TBSA, which converts TBSA to a tagging task, and aims to solve the two subtasks (i.e. target detection and sentiment classification) in a unified manning by predicting unified tags. For example, the tag ``B-POS'' indicates the beginning of a target with positive sentiment. So after annotation, the above example becomes \emph{``[B-POS USB3] [E-POS Peripherals] are noticeably less expensive than the [B-NEG ThunderBolt] [E-NEG ones]"}.

\section{Proposed Method}
\label{sec:method}

We propose a novel data augmentation method for sequence tagging tasks. We first linearize labeled sentences, and then train a language model to learn the distribution of words and tags from the linearized sequences for generating synthetic training data. A conditional generation technique is also proposed to exploit unlabeled data and knowledge bases when they are available.

\subsection{Labeled Sentence Linearization}
\label{ssec:sentence_linearization}
We first perform sentence linearization to convert labeled sentences into linear sequences, so that language models can be used to learn the distribution of words and tags in gold data. As shown in Figure~\ref{fig:linearization}, tags are inserted before the corresponding words during linearization and thus treated as modifiers of these words. For the tasks with frequent \emph{O} tags,  e.g., NER and E2E-TBSA \citep{li2019unified}, we remove such tags from the linearized sequences. Similarly, we can also insert tags after the corresponding words.

After sentence linearization, we add special tokens \emph{[BOS]} and \emph{[EOS]} to the beginning and the end of each sentence,  respectively.
These special tokens are used to facilitate model training and data generation by marking the sentence boundaries.

\subsection{Language Modeling and Data Generation}
\label{subsec:lm_and_gen}
After linearizing labeled sentences, language models can be used to learn the distribution of words and tags.  More specifically, we use a one-layer LSTM recurrent neural network language model (RNNLM) in our method, which is similar to the model proposed by \citet{sundermeyer2012lstm}. The architecture of our RNNLM is shown in Figure~\ref{fig:rnnlm}.

\begin{figure}[t!]
\centering
\includegraphics[scale=0.38]{image/rnnlm.png}
\caption{Language model architecture with LSTM.}
\label{fig:rnnlm}
\end{figure}

\paragraph{Language Modeling}

We train RNNLM by maximizing the probability for next token prediction. Given a sentence, we first feed the sequence of tokens $(w_1, w_2, \dots, w_N)$ into the embedding layer to lookup the corresponding embeddings $({\bm{x}}_1, {\bm{x}}_2, \dots, {\bm{x}}_N)$, where $N$ is the sequence length. A dropout layer is applied to each token embedding ${\bm{x}}_t$ to generate ${\bm{d}}_t = \mathrm{dropout}({\bm{x}}_t)$. Then we feed $({\bm{d}}_1, {\bm{d}}_2, \dots, {\bm{d}}_N)$ into LSTM to produce hidden state ${\bm{h}}_t = \mathrm{LSTM}({\bm{d}}_t,{\bm{h}}_{t-1})$ at each position $t$. Another dropout layer is applied to hidden states to compute ${\bm{d}}^\prime_t = \mathrm{dropout}({\bm{h}}_t)$.

\begin{figure*}[t!]
\centering
\includegraphics[scale=0.45]{image/conditional_generation.png}
\caption{An example of conditional generation. The first sequence is from gold NER data. The second is from unlabeled data, so no labels. The third is labeled by knowledge base matching, where \emph{Asakusa} cannot be labeled due to incomplete knowledge coverage.}
\label{fig:conditional_generation}
\end{figure*}

Finally, a linear+softmax layer is used to predict the next token in the sequence. Assuming the index of token $w_t$ in the vocabulary is $i^*$, we have the training objective in Eq.~\ref{eq:train_objective}:
\begin{eqnarray}
    {\bm{s}}_{t-1} = {\bm{M}}^\mathsf{T}{\bm{d}}^\prime_{t-1} \\

    p_{\theta}(w_t|w_{<t}) = \frac{\exp({\bm{s}}_{t-1,i^*})}{\sum_{i=1}^V \exp({\bm{s}}_{t-1,i})}
    \label{eq:word_prediction} \\
    p(w_1, w_2, \dots, w_N) = \prod_{t=1}^N p_{\theta}(w_t|w_{<t})
    \label{eq:train_objective}
\end{eqnarray}
where $V$ is the size of vocabulary, ${\bm{M}} \in \mathbb{R}^{r \times V}$ is a learnable weight matrix with $r$ being the dimension of LSTM hidden states, and ${\bm{s}}_{t-1,i}$ is the $i$-th element of ${\bm{s}}_{t-1}$.

\paragraph{Generation}
After training the RNNLM, we can use it to generate synthetic training data for tagging tasks. During generation, only the \emph{[BOS]} token is fed into RNNLM, and the following tokens are sampled based on the probabilities computed by Eq.~\ref{eq:word_prediction}. Given \emph{[BOS]}, the sentence is generated autoregressively one at a time, where the token generated in the previous step is taken as input to generate the next one.

As shown in Eq.~\ref{eq:word_prediction}, RNNLM is more likely to pick the tokens with high probabilities during sampling in the generation process. Because of the randomness added by sampling, RNNLM can choose similar alternatives given the same context. Assume we insert tags before the corresponding words (during sentence linearization) to train the RNNLM, when predicting the next token given \emph{``I have booked a flight to"}, the probability of \emph{``S-LOC"} is much higher than the other choices, since the RNNLM has seen many similar examples in training data, such as \emph{``a train to S-LOC"}, \emph{``a trip to S-LOC"} and so on. Then we predict the following word given \emph{``I have booked a flight to S-LOC"}. In the training data, all \emph{``S-LOC"} are followed by location words, so \emph{``London"}, \emph{``Paris"}, \emph{``Tokyo"}, etc., are all possible choices, and their probabilities are very close. Due to the added randomness, the model can choose any one of them. Tokens are predicted in a similar way when we insert tags after the corresponding words, except that words are predicted before the tags.

\subsection{Post-Processing}
\label{sec:post_process}
The generated sequences are in the linearized format, so they need to be converted to the same format as the gold data. We also introduce several straightforward rules to clean the generated data: 1) Delete sentences with no tags; 2) Delete sentences where all words are \emph{[unk]}\footnote{To reduce the size of vocabulary when training language model, the words that only appear once in training data are replaced with the unknown token \emph{[unk]}.}; 3) Delete sentences with incorrect tag prefix orders (e.g., having \emph{E-LOC} before \emph{B-LOC} in NER data); 4) Delete sentences that contain same sequences of words but different tags.

\subsection{Conditional Generation}
\label{subsec:conditional_generation}
We propose a conditional generation method to allow the language model to utilize unlabeled data or knowledge bases when they are available in some low-resource scenarios. For example, it could be expensive to annotate a large amount of e-commerce product titles for NER, but much easier to obtain a knowledge base (i.e., dictionary) of product attributes and unlabeled data. We prepend one of these \textbf{condition tags} $\{ \emph{[labeled]}, \emph{[unlabeled]}, \emph{[KB]} \}$ at the beginning of each sequence to mark their origin, where \emph{KB} means the sequence is labeled by matching a knowledge base against the unlabeled data. See Figure~\ref{fig:conditional_generation} for an example. This allows the language model to learn the shared information among these sequences while being aware of the different origins. When generating synthetic data, each word is conditioned on the given condition tag \emph{[labeled]}, denoted as $c$ (conditioning class). After we feed it into the language model, all the LSTM hidden states ${\bm{h}}$ in the following generation steps contain information of $c$. In addition, ${\bm{h}}$ also encodes information of all of the other previous tokens in the sequence. When predicting the next token conditioned on ${\bm{h}}_{t-1}$, the probability $p_{\theta}(w_t|w_{<t})$ in Eq.~\ref{eq:word_prediction} becomes $p_{\theta}(w_t|w_{<t},c)$.\footnote{Condition tag $c$ is also in $w_{<t}$, since it is a special token added to the beginning of each sentence. We write it explicitly to emphasize the conditional effect.} A similar approach
is used in CTRL \citep{keskar2019ctrl} to control style, task-specific behavior, etc., during text generation.

\section{Experiments} \label{sec:exp}

In this section, we present our experiments in both supervised and semi-supervised settings. In the supervised settings, only gold data are used for augmentation. In the semi-supervised settings, we also leverage unlabeled data and knowledge bases.

\subsection{Basic Models}
\paragraph{Language Model}
We use the language model described in Section~\ref{subsec:lm_and_gen} for synthetic data generation. We modified the decoder of the LSTM-LM model in \citet{kruengkrai2019better} to implement this language model. We set LSTM hidden state size to 512 and embedding size to 300. We use dropout rate 0.5 for the two dropout layers. All language model are trained using Stochastic gradient descent (SGD) with initial learning rate 1 and batch size 32. Learning rate will be decayed by 0.5 in the next epoch if the perplexity on dev set does not improve. We set the maximum number of epochs to 30 and stop training early if the perplexity on dev set does not improve in 3 consecutive epochs. During synthetic data generation, we use the average length of gold sentences in the training set as our maximum sentence length.

\paragraph{Sequence Tagging Model}
We implement a BiLSTM-CRF model \citep{lample2016neural} with the Flair framework \citep{akbik2019flair} to evaluate our data augmentation method on NER and POS tasks.\footnote{The baseline model provided in the original paper is used for evaluating end to end target based sentiment analysis task.} We use a single-layer BiLSTM with hidden state size 512. Dropout layers are applied before and after the BiLSTM layer with dropout rate 0.5. All sequence tagging models are trained using Adam \citep{kingma2014adam} with initial learning rate 1e-3 and batch size 32. Learning rate is decayed by 0.5 if the performance on dev set does not improve in 3 consecutive epochs. We stop training when the learning rate drops below 1e-5 or number of epochs reaches 100. We use the pre-trained 300-dimensional fastText word embeddings \citep{bojanowski2017enriching} for all languages.

We employ relatively simple basic models because: 1) They help to avoid the possible overfitting problems due to the small data size under the low resource setting; 2) They allow more faithful understanding on the effects of the proposed data augmentation method.

\label{sec:experiment}
\subsection{Supervised Experiments}
To verify the effectiveness of our data augmentation method in the supervised settings, we evaluate it on three different tagging tasks, including NER, POS and E2E-TBSA. Most of the prior works rely on additional information, so we use random deletion (\textbf{rd}) \citep{wei2019eda} as our baseline, where 5\% of the words\footnote{For NER and E2E-TBSA, the whole entity is deleted if a selected word appears within an entity span.} and the corresponding tags in training data are randomly deleted. See Table~\ref{table:supervised_approaches} for the notations of the methods used in our supervised experiments.

\begin{table}[!ht]
\centering
\scalebox{0.8}{\begin{tabular}{@{}l@{~}@{~}p{8cm}@{}}
\toprule
\textbf{Method} & \textbf{Description}\\
\midrule
\textbf{gold} & Only use the gold data.  \\
\textbf{gen} & Our method. Generate synthetic data with the language models, and oversample gold data. \\
\textbf{rd} & Baseline method. Generate synthetic data by random deletion, and oversample gold data with the same ratio as \textbf{gen}. \\
\textbf{rd\textsuperscript{*}} & Baseline method. Similar to \textbf{rd}, except that gold and synthetic data are equally sampled.  \\
\bottomrule
\end{tabular}}
\caption{Data sources for the supervised setting.}
\label{table:supervised_approaches}
\end{table}

\subsubsection{Named Entity Recognition}

\paragraph{Dataset}
We evaluate our proposed methods on the CoNLL2002/2003 NER data \citep{tjong-kim-sang-2002-introduction,tjong-kim-sang-de-meulder-2003-introduction}, with four languages: English, German, Dutch and Spanish. Besides, we evaluate our methods on Thai and Vietnamese NER data, which are product titles obtained from major e-commerce websites in Southeast Asian countries and annotated with 11 product attribute NER tags, including \emph{PRODUCT}, \emph{BRAND}, \emph{CONSUMER\_GROUP}, \emph{MATERIAL}, \emph{PATTERN}, \emph{COLOR}, \emph{FABRIC}, \emph{OCCASION}, \emph{ORIGIN}, \emph{SEASON} and \emph{STYLE}. See Appendix for the statistics of the Thai and Vietnamese NER data used in our experiments.

\paragraph{Experimental Settings}
In addition to evaluating our method on the full training data, we also randomly sample 1k, 2k, 4k, 6k and 8k sentences for each language to verify its robustness in low-resource settings. We use the original development and test data. The language models are trained on the above sampled sentences, and then used to generate synthetic training data following the steps described in Section~\ref{sec:method}.

For each new batch of 1k sentences generated by the trained language model, we measure the percentage of new 1-gram tokens that appears in previous batches. Once the percentage exceeds 99\%, we will stop data generation. Then we post-process (Section \ref{sec:post_process}) the generated data, and add them to gold training data for tagging model training. For \textbf{rd} and \textbf{gen}, we oversample gold data by repeating them 4 times (shuffled) in the training set. For random deletion, we also report the result when gold and synthetic data are equally sampled, denoted by \textbf{rd\textsuperscript{*}}. Additional details on hyperparameter selection for the oversampling ratios can be found in Appendix A.3. Following \citet{lample2016neural}, the IOBES tagging scheme is used when training the language models and sequence tagging models described above.

\paragraph{Results and Analysis}
We report the average results of 3 runs in Table~\ref{table:conll_ner}. Our method shows consistent performance improvement for all languages. Especially for the smaller sampled sets, our method demonstrates more significant performance improvement. In particular, the proposed method achieved average 1.93 and 1.38 point improvement compared with the baseline methods in the 1k and 2k settings, respectively.

\begin{table}[t!]
\centering
\scalebox{0.75}{\begin{tabular}{@{~}cl@{~}ccccccc@{~}}
\toprule
\textbf{Lang.} & \textbf{Method} & \textbf{1k} & \textbf{2k} & \textbf{4k} & \textbf{6k} & \textbf{8k} & \textbf{all} \\
\midrule
\multirow{4}{*}{\textbf{en}} & gold & 58.06 & 67.85 & 74.55 & 77.16 & 80.30 & 83.04 \\
& ~+rd\textsuperscript{*} & 59.42 & 67.23 & 74.51 & 77.39  & 80.31  & 83.39          \\
& ~+rd & 58.97 & 67.81 & 74.77 & 77.35 & 80.59  & 83.25          \\
& ~+gen & \textbf{61.15} & \textbf{70.61} & \textbf{76.82} & \textbf{79.18} & \textbf{81.02} & \textbf{83.74} \\
\midrule
\multirow{4}{*}{\textbf{de\tablefootnote{German has large out-of-vocabulary rate when using static fastText embedding, which leads to lower F1 score compared with other languages.}}} & gold & 29.71  & \textbf{41.07} & 49.55 & 53.30 & 56.17 & 61.10 \\
& ~+rd\textsuperscript{*} & 29.89  & 40.29  & 49.27  & 52.33  & 55.70 & 60.69 \\
& ~+rd & 30.83 & 40.36   & 49.24   & 53.54  & 55.60   & 60.55          \\
& ~+gen &\textbf{31.83} & 40.92 & \textbf{49.79} & \bf{53.63} & \textbf{56.94} & \textbf{62.44} \\
\midrule
\multirow{4}{*}{\textbf{es}} & gold & 58.14  & 67.42  & 74.21 & 77.44 & 78.90 & 79.27          \\
& ~+rd\textsuperscript{*} & 58.22  & 66.98  & 75.08   & 77.64   & 79.11   & 80.01          \\
& ~+rd & 59.67  & 68.53  & 75.21 & 77.79  & 79.12    & 80.26          \\
& ~+gen & \textbf{61.76} & \textbf{68.62} & \textbf{76.15} & \textbf{78.20} & \textbf{79.83} & \textbf{80.73}  \\
\midrule
\multirow{4}{*}{\textbf{nl}} & gold & 37.04 & 48.61 & 57.78 & 61.08 & 64.59 & 70.89   \\
& ~+rd\textsuperscript{*} & 35.10    & 46.45   & 56.83  & 60.49  & 63.09   & 69.42          \\
& ~+rd & \textbf{39.39} & 48.44   & 59.38  & 61.48 & 64.44   & 70.36          \\
& ~+gen & 38.87  & \textbf{50.41} & \textbf{59.90} & \textbf{63.19} & \textbf{65.82} & \textbf{72.71}  \\
\midrule
\multirow{4}{*}{\textbf{vi}} & gold & 55.98  & 62.42  & 69.01 & 70.75  & 72.12 & 76.14       \\
& ~+rd\textsuperscript{*} & 55.67 & 63.57 & 68.47 & 70.87 & 72.08 & 76.43 \\
& ~+rd & 56.24 & 63.08 & 68.63 & 71.15 & 72.22 & 76.83 \\
& ~+gen & \textbf{60.01} & \textbf{65.43} & \textbf{70.36} & \textbf{72.55} & \textbf{74.11} & \textbf{77.39} \\
\midrule
\multirow{4}{*}{\textbf{th}} & gold & 49.88 & 55.79 & 61.75 & 63.10 & 64.94 & 67.71   \\
& ~+rd\textsuperscript{*} & 50.46 & 56.98 & 62.12 & 64.19 & 66.47 & 67.81 \\
& ~+rd & 50.52 & 57.42 & 61.51 & 64.59 & 66.07 & 67.97 \\
& ~+gen & \textbf{54.02} & \textbf{59.36} & \textbf{63.94} & \textbf{66.21} & \textbf{68.05} & \textbf{69.86} \\
\bottomrule
\end{tabular}}
\caption{Named entity recognition micro F1.}
\label{table:conll_ner}
\end{table}

\paragraph{Tag-Word vs. Word-Tag}
As discussed in Section~\ref{ssec:sentence_linearization}, there are two ways to perform sentence linearization: 1) insert tags before the corresponding words (Tag-Word); 2) insert tags after the corresponding words (Word-Tag). Keeping all of the other settings same, we find Tag-Word outperforms Word-Tag in NER tasks (as shown in Appendix A.2). One possible reason is that, Tag-Word is more consistent with the Modifier-Noun pattern, which appears more often in the training data during language modeling. Therefore, we use Tag-Word for all the NER experiments.

\subsubsection{Part of Speech Tagging}
\paragraph{Dataset}
We use the POS data from Universal Dependencies treebanks\footnote{https://universaldependencies.org/} for evaluation on this task. We evaluate on five languages, including English, Spanish, Czech, Romanian and Japanese. Each language has multiple corpora in the Universal Dependencies treebanks, so we merge the corpora to build one dataset for three languages: English, Spanish and Czech. For English, we merge \emph{GUM}, \emph{ParTUT}, \emph{PUD} and \emph{Lines}. For Spanish, we merge \emph{AnCora} and \emph{GSD}. For Czech, we merge \emph{PDT}, \emph{FicTree}, \emph{CLTT} and \emph{CAC}. We have also evaluated our model on Japanese (\emph{GSD}) and Romanian (\emph{RRT}), which are either spoken by a much smaller population or from a different language family.

\paragraph{Settings and Results}
We follow similar experimental settings as NER task. The same language model and BiLSTM-CRF sequence tagging model are used for synthetic data generation and POS tagging respectively. Different from NER, Word-Tag shows slightly better performance in POS tasks (refer to Appendix A.2).
We present the average Word-Tag results of 3 runs in Table~\ref{table:ud_pos}. Our method demonstrates consistent performance improvement for all languages. Similar to NER, our method demonstrates more significant performance improvement for smaller sampled sets on POS tagging. In particular, the proposed method achieved average 0.56, 0.60 and 0.46 point improvement compared with the baseline methods in the 1k, 2k and 4k settings, respectively.

\begin{table}[t!]
\centering
\scalebox{0.75}{\begin{tabular}{@{~}cl@{~}cccccc@{~}}
\toprule
\textbf{Lang.} & \textbf{Method} & \textbf{1k} & \textbf{2k} & \textbf{4k} & \textbf{6k} & \textbf{8k} & \textbf{Full} \\
\midrule
\multirow{4}{*}{\textbf{en}} & gold & 79.18 & 82.17 & 85.83 & 88.62 & 90.21  & 93.00\\
& ~+rd\textsuperscript{*} & 79.28 & 82.42 & 85.82 & 88.55 & 90.07 & 92.89 \\
& ~+rd & 79.38 & 82.50 & 86.08 & 88.80 & 90.15 & 92.96 \\
& ~+gen & \bf{79.76} & \bf{82.90} & \bf{86.31} & \bf{88.99} & \bf{90.56} & \bf{93.29} \\
\midrule
\multirow{4}{*}{\textbf{es}} & gold & 88.28 & 90.79 & 92.82 & 93.80 & 94.43  & 96.40 \\
& ~+rd\textsuperscript{*} & 88.25 & 90.94 & 92.84 & 93.76 & 94.48 & 96.41\\
& ~+rd & 88.17 & 90.78 & 92.79 & 93.67 & 94.28 & \bf{96.45}\\
& ~+gen & \bf{88.77} & \bf{91.04} & \bf{93.12} &	\bf{93.93} & \bf{94.64} & \bf{96.45}\\
\midrule
\multirow{3}{*}{\textbf{cz}} & gold & 80.10 & 84.46 & 88.88 & 90.67 & 92.03 & 97.52\\
& ~+rd\textsuperscript{*} & 79.83 & 84.29 & 88.64 & 90.43 & 91.95 & 97.57 \\
& ~+rd & 80.11 & 84.50 & 88.99 & 90.66 & 91.86 & 97.60 \\
& ~+gen & \bf{80.65} & \bf{85.17} & \bf{89.58} & \bf{91.22} & \bf{92.49} & \bf{97.63} \\
\midrule
\multirow{3}{*}{\textbf{ro\tablefootnote{UD-RRT full train set has 8k sentences for Romanian.}}} & gold & 86.69 & 89.57 & 92.73 & 93.84 & 94.54 & 94.54 \\
& ~+rd\textsuperscript{*} & 86.42 &	89.58 &	92.50 &	93.89 &	94.64 & 94.64 \\
& ~+rd & 86.62 & 89.46 & 92.55 & 93.84 & 94.73 & 94.73  \\
& ~+gen & \bf{87.29} &	\bf{90.66} & \bf{93.44} & \bf{94.61} & \bf{95.17} & \bf{95.17}\\
\midrule
\multirow{3}{*}{\textbf{ja\tablefootnote{UD-GSD full train set has 7k sentences for Japanese.}}} & gold & 90.19 & 91.44 & 93.59 & 94.41 & - & 95.08
\\
& ~+rd\textsuperscript{*} & 90.00 &	91.41 &	93.66 &	94.62 & - &	94.93 \\
& ~+rd & 89.53 & 91.76 & 93.62 & 94.59 & - 	& 95.18 \\
& ~+gen & \bf{91.00} & \bf{92.51} &	\bf{94.12} & \bf{95.21} & - & \bf{95.45} \\
\bottomrule
\end{tabular}}
\caption{POS tagging accuracy.}
\label{table:ud_pos}
\end{table}

\subsubsection{Target Based Sentiment Analysis}
\paragraph{Dataset}
We use the laptop and restaurant review datasets processed by \citet{li2019unified} for evaluation on E2E-TBSA, which was initially obtained from SemEval ABSA challenges \citep{pontiki-etal-2014-semeval,pontiki-etal-2015-semeval,pontiki-etal-2016-semeval}. We merge these two review datasets, regard $10\%$ randomly held-out training data as the dev set, and randomly sample smaller sets from the remaining training data for low-resource settings. The original test sets are merged as our test set.

\paragraph{Settings and Results}
We follow similar experimental settings as NER and POS tasks, except that the same sequence tagging model released by \citet{li2019unified} is used for evaluation. Here  Tag-Word shows better results (refer to Appendix A.2), plausibly it is because the unified tags (e.g. \emph{B-POS}, and \emph{B-NEG}) are similar to noun modifiers and Tag-Word is more consistent with the Modifier-Noun pattern. We present the average Tag-Word results of 3 runs in Table~\ref{table:e2etbsa}. Our method demonstrates performance improvement for 4k and above. Compared with the NER and POS datasets, the E2E-TBSA dataset has much fewer labels, so the results are less stable.

\begin{table}[th!]
\centering
\scalebox{0.78}{\begin{tabular}{lcccc}
\toprule
\textbf{Method} & \textbf{2k} & \textbf{4k} & \textbf{all(6k)} \\
\midrule
gold & 56.31 & 60.43 & 63.18  \\
~+rd\textsuperscript{*} & \bf{57.92} & 61.75 & 63.66 \\
~+gen & 57.07 & \bf{62.66} & \bf{65.86}  \\
\bottomrule
\end{tabular}}
\caption{ E2E-TBSA micro F1.}
\label{table:e2etbsa}
\end{table}

\subsection{Semi-supervised Experiments}
In this section we evaluate the effectiveness of our method in two semi-supervised settings: 1) only unlabeled data are available; 2) both unlabeled data and knowledge base are available. See Table~\ref{table:semisupervised_approaches} for the notations of the methods used in our semi-supervised experiments.

\begin{table}[htb!]
\centering
\scalebox{0.8}{\begin{tabular}{@{}l@{~}@{~}p{8cm}@{}}
\toprule
\textbf{Method} & \textbf{Description}\\
\midrule
\textbf{gold} & Supervised method. Only use the gold data.  \\
\textbf{wt} & Baseline method. Annotate unlabeled data with a weak tagger (i.e. a tagging model trained on the gold data).  \\
\textbf{gen\textsubscript{ud}} & Our method. Generate synthetic data with LM, where LM is trained on gold data and unlabeled data. \\
\textbf{kb} & Baseline method. Annotate unlabeled data with knowledge base.  \\
\textbf{gen\textsubscript{kb}} & Our method. Generate synthetic data with LM, where LM is trained on gold data and knowledge base annotated data. \\
\bottomrule
\end{tabular}}
\caption{Data sources for the semi-supervised setting.}
\label{table:semisupervised_approaches}
\end{table}

\subsubsection{Only Using Unlabeled Data}
\label{sssec:semi_ner_unlabeled}

\paragraph{Dataset}
We use CoNLL2003 English NER data \citep{tjong-kim-sang-de-meulder-2003-introduction} for evaluation. In addition to the gold NER training data, we utilize unlabeled data for semi-supervised training. The Stanford CoreNLP tokenizer \citep{manning-EtAl:2014:P14-5} is used to tokenize Wikipedia sentences.

\paragraph{Experimental Settings}
Similar to the above experiments, we use 1k, 2k, 4k, 6k and 8k sentences randomly sampled from NER gold data as well as the full dataset to evaluate our method. For fair comparison, we only use the same set of 10k sentences randomly sampled from Wikipedia dump in both of our and baseline methods. Let $D_{gold}$ and $D_{unlabeled}$ be the sampled gold NER data and the Wikipedia data, respectively.

In our method, $D_{gold}$ and $D_{unlabeled}$ are concatenated to train language models, following the steps described in Section~\ref{subsec:conditional_generation}.
Then we use the language models to generate synthetic data, from which 20k randomly sampled sentences are combined with $D_{gold}$ to train NER models. We use \textbf{gen\textsubscript{ud}} to denote this data generation method. The method that employs weak taggers to annotate $D_{unlabeled}$ is used as our baseline, denoted by \textbf{wt}. The weak taggers in this experiment are the NER models trained on $D_{gold}$. We use the same NER model (BiLSTM-CRF) and hyperparameters to evaluate our and baseline methods. When training the language model, we equally sample sentences from $D_{gold}$ and $D_{unlabeled}$. When training the NER models, we oversample the gold data by repeating $D_{gold}$ 4 times to create a shuffled training file.

\paragraph{Results and Analysis}
We report F1 of \textbf{wt} and \textbf{gen\textsubscript{ud}} (average of 3 runs) in Table~\ref{table:semi_ner_unlabeled}. Our method outperforms the baseline method \textbf{wt} on all settings. Moreover, it  would be a promising direction to further explore our model's capability by using a larger amount of unlabeled data. It is not very convenient to utilize a large amount of unlabeled data in the baseline method \textbf{wt}, since this will directly increase the amount of augmented data. As a result, some of augmented data may not be utilized before sequence tagging models converge. However, our method \textbf{gen\textsubscript{ud}} can conveniently utilize a large amount of unlabeled data to train the language models, thereby improving data augmentation quality directly. When the amount of unlabeled data is much larger, we can pretrain language models with the unlabeled data, and then finetune them with labeled data.

\begin{table}[t!]
\centering
\scalebox{0.82}{\begin{tabular}{lcccccc}
\toprule
\textbf{Method} & \textbf{1k} & \textbf{2k} & \textbf{4k} & \textbf{6k} & \textbf{8k} & \textbf{all} \\
\midrule
gold & 58.06 & 67.85 & 74.55 & 77.16 & 80.30 & 83.04 \\
\midrule
\rowcolor{lightyellow}  ~+wt & 65.12 & 72.43 & 77.90 & 79.41 & 81.36 & 84.00 \\
\rowcolor{lightyellow}  ~+gen\textsubscript{ud} & \bf{66.19} & \bf{73.00} & \bf{78.08} & \bf{79.75} & \bf{81.98} & \bf{84.33} \\
\midrule
\rowcolor{lavender}~+kb & \bf{67.36} & 72.86 & 77.15 & 79.33 & 81.91 & 83.69 \\
\rowcolor{lavender}~+gen\textsubscript{kb} & 66.67 & \bf{73.54} & \bf{78.32} & \bf{79.98} & \bf{81.93} & \bf{84.03} \\
\bottomrule
\end{tabular}}
\caption{Semi-supervised NER F1.}
\label{table:semi_ner_unlabeled}
\end{table}

\begin{comment}
\begin{table}[hbt!]
\centering
\scalebox{0.78}{\begin{tabular}{lcccccc}
\toprule
\textbf{Method} & \textbf{1k} & \textbf{2k} & \textbf{4k} & \textbf{6k} & \textbf{8k} & \textbf{all} \\
\midrule
gold & 58.06 & 67.85 & 74.55 & 77.16 & 80.30 & 83.04 \\
~+wt & 65.12 & 72.43 & 77.90 & 79.41 & 81.36 & 84.00 \\
~+gen\textsubscript{ud} & \bf{66.19} & \bf{73.00} & \bf{78.08} & \bf{79.75} & \bf{81.98} & \bf{84.33} \\

\bottomrule
\end{tabular}}
\caption{Semi-supervised NER F1 (unlabeled data)}
\label{table:semi_ner_unlabeled}
\end{table}
\end{comment}

\subsubsection{Using Unlabeled Data and Knowledge Base}
\paragraph{Dataset}
In addition to the gold training data and unlabeled sentences used in Section~\ref{sssec:semi_ner_unlabeled}, we also try to leverage knowledge base for further performance improvement in this experiment. We build the knowledge base by extracting entities (case sensitive and appearing at least twice) and the corresponding tags from the full gold training data. Besides, we add more \emph{LOC} entities to this knowledge base by including the cities and countries extracted from geonames\footnote{https://datahub.io/core/world-cities}.

\paragraph{Experimental Settings}
We randomly sample the gold NER data and the Wikipedia data in the same way as Section~\ref{sssec:semi_ner_unlabeled}, where the sampled sentences are denoted by $D_{gold}$ and $D_{unlabeled}$, respectively. Our knowledge base is used to annotate $D_{unlabeled}$ by finding longest forward matches (from left to right) in each sentence. We denote this annotation method by \textbf{kb} and the annotated data by $D_{kb}$.

In our method, $D_{gold}$ and $D_{kb}$ are concatenated to train language models, following the steps described in Section~\ref{subsec:conditional_generation}.
Then we use the language models to generate synthetic data, from which 20k randomly sampled sentences are combined with $D_{gold}$ to train the NER models.
We use \textbf{gen\textsubscript{kb}} to denote this data generation method and compare it with the baseline method \textbf{kb}. Similar to the above experiments, we oversample $D_{gold}$ when training the language and NER models.

\paragraph{Results and Analysis}
We present F1 of \textbf{kb} and \textbf{gen\textsubscript{kb}} (average of 3 runs) in Table~\ref{table:semi_ner_unlabeled}. The baseline method \textbf{kb} exhibits very strong performance when the size of $D_{gold}$ is small, since we use a large knowledge base of countries and cities for annotation, and location names are less ambiguous compared with the other types of entities. However, our method still outperforms \textbf{kb} when the size of $D_{gold}$ is larger than 2k, which shows that our method is more robust to the noises in $D_{kb}$ when a slightly larger amount of gold data are available.

\begin{comment}
\begin{table}[h!]
\centering
\scalebox{0.78}{\begin{tabular}{lcccccc}
\toprule
\textbf{Method} & \textbf{1k} & \textbf{2k} & \textbf{4k} & \textbf{6k} & \textbf{8k} & \textbf{all} \\
\midrule
gold & 58.06 & 67.85 & 74.55 & 77.16 & 80.30 & 83.04 \\
~+kb & \bf{67.36} & 72.86 & 77.15 & 79.33 & 81.91 & 83.69 \\
~+gen\textsubscript{kb} & 66.67 & \bf{73.54} & \bf{78.32} & \bf{79.98} & \bf{81.93} & \bf{84.03} \\
\bottomrule
\end{tabular}}
\caption{Semi-supervised NER F1 (unlabeled data and knowledge base)}
\label{table:semi_ner_kb}
\end{table}
\end{comment}

\section{A Closer Look at Synthetic Data}
\label{sec:gen_data}

In this section, we explore in more details why the synthetic data generated by our method can help improve sequence tagging performance. Through a closer look at the generated data, we have several interesting findings.

\paragraph{More Diversity}
The generated synthetic data introduces more diversity to help reduce overfitting. As the example shown in Figure~\ref{fig:analysis_diversity}, the name ``\emph{Sandrine}" in the gold training data always pairs up with ``\emph{Testud}" in different sentences. However, in the generated data, we can see new names have been generated like ``\emph{Sandrine Nixon}", ``\emph{Sandrine Okuda}" and ``\emph{Sandrine Neuumann}". Meanwhile, the locations in  the sentences have been replaced with new countries like ``\emph{Sweden}", ``\emph{Egypt}" and ``\emph{Australia}". With these synthetic data, the model can focus on learning the pattern of contexts that entities appear in, instead of simply memorizing ``\emph{Sandrine Testud}" as a person name and ``\emph{France}" as a location.

To quantitatively measure the diversity brought in by our method and its impact, we have done some statistical analysis of the contextualized entities (CEs) in the generated data of the supervised English NER. A CE refers to the combination of an entity and its 1-gram contexts. For example, in the sentence ``\emph{The [B-ORG European] [E-ORG Commission] said ...}", the entity is ``\emph{European Commission}" and its CE is ``\emph{The European Commission said}".
As the shown in Figure~\ref{fig:contextualized_entity}, we calculate the number of unique CEs in the gold training data, the number of new unique CEs in the generated data, and the ratio of the two numbers. We also plot here the F1 improvement of our method (i.e. gold+gen) over only using the gold data, as given in Table \ref{table:conll_ner}.
We can see that our method generates a large number of new CEs, and such diversity strengthens the robustness of the trained model.
When the ratio is higher, the F1 improvement is more significant, which also shows that our method does help ease the low-resource problem by generating rich new entities and contexts.
Refer to the Appendix A.5 for statistics of unique entities (without context), which shows the same conclusion.

\begin{figure}[t!]
\centering
\includegraphics[scale=0.4]{image/generated_data_new.png}
\caption{An illustration of diversity of generated data. The name ``\emph{Sandrine}" in the gold training data always pairs up with ``\emph{Testud}" in sentences.}
\label{fig:analysis_diversity}
\end{figure}

\paragraph{Efficient Utilization of Unlabeled Data}
When unlabeled data are available, our method is flexible to utilize them for semi-supervised training. We find many interesting examples in the synthetic data that show our method can effectively use the unlabeled data to extract useful information. In an example generated by our method ``\emph{... the [B-ORG Bank] [I-ORG of] [E-ORG Alabama] ...}", the word ``\emph{Alabama}” has never appeared in gold NER training data. However, our language model learned that ``\emph{Alabama}” (from unlabeled data) is very similar to the other location words that appear in both gold training data and unlabelled data. So when generating the synthetic data, the language model can use this word in a similar context, or even create new entities (``\emph{Bank of Alabama}" in this example has never appeared in gold or unlabeled data).

\begin{figure}[t!]
\centering
\includegraphics[scale=0.52]{image/Diversity.png}
\caption{Statistics of unique contextualized entities.}
\label{fig:contextualized_entity}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we show that language models can be used to generate high quality synthetic data for sequence tagging tasks. The generated data introduce more diversity to reduce overfitting, since they are generated from scratch instead of modifying the gold training data. Our proposed method demonstrates promising performance improvements on various tagging tasks, especially in the low-resource settings. Besides, experiments demonstrate that our method can also effectively utilize unlabeled data and knowledge base for semi-supervised training.

\section*{Acknowledgements}
\label{sec:acknowledgements}

This research is partly supported by the Alibaba-NTU Singapore Joint Research Institute, Nanyang Technological University. Linlin Liu would like to thank the support from Interdisciplinary Graduate School, Nanyang Technological University.

\bibliographystyle{acl_natbib}


\section{ Introduction }
Let $\Omega_g=\Omega_2\times (0,g)$ where $\Omega_2$ is a bounded domain in $\mathbb R^2$ and $g$ is some scalar nonnegative function. We introduce the following time-fractional $g$-B\'enard equations of the following form:
  \begin{align*}
  \partial^{\alpha}_t u+(u\cdot\nabla)\,u-\nu  \Delta u+\nabla p&= \xi \theta +f_1(t)  \\[0.1cm]
  \nabla\cdot gu &= 0\\[0.1cm]
  \partial^{\alpha}_t \theta+(u\cdot\nabla)\,\theta-\kappa\Delta\theta &= f_2(t)
  \end{align*}
where $u$ is the fluid velocity, $p$ is the pressure, $\theta$ is the temperature, $f_1$ is the external force function, $f_2$ is the heat source function, $\xi\in\mathbb{ R}^3$ is a constant vector, $\nu$ the kinematic viscosity and $\kappa$ the thermal diffusivity are positive constants. The derivative of order $\alpha$ is considered in the Caputo sense. The time-fractional $g$-B\'enard problem consists in a system that couples time-fractional Navier-Stokes equations and time-fractional advection-diffusion heat equation in order to model a memory-dependent convection in a fluid considered in a fractal media.

The introduction of time-fractional derivative in fluid dynamics goes back to Lions in \cite{L59} but for order less than $\frac{1}{4}$ provided the space dimension is not further than $4$.  In recent works of Zhou and Peng \cite{ZP17, ZP172}, the question of weak solutions and optimal control problem of time-fractional Navier-Stokes equations in fractal media was considered. Numerical results regarding such problems was treated firstly in \cite{KKAR14} and constitute an emerging field of research. More recently, time-fractional $g$-Navier-Stokes problem was introduced and results regarding existence, uniqueness of solutions and optimal control are proved\cite{BA20}.

The theory of $g$-Navier-Stokes equations started with the works of Hale and Raugel \cite{HR92a, HR92b}, Raugel and Sell \cite{RS93} who studied 3d nonlinear equations and Navier-Stokes equations in thin domains. J. Roh\cite{R01}, a student of Sell, generalized the previous works to thin domains of the form $ \Omega_g=\Omega\times(0,g)$, where $g$ is some smooth scalar function. The derived equations are called the $g-$Navier-Stokes equations. This theory have interested many researchers in recent years, see  \cite{BR04,AQ12a,JHW11} and references therein.

On the other hand, heat conduction based on the classical Fourier law, which relates the heat flux vector and the temperature gradient, has shown its limits. The time-fractional heat conduction model can be seen as a good alternative\cite{MLP01,GIL00,SZ97}, and references therein. Boussinesq(or B\'enard) model is a combination of the heat conduction model and Navier-Stokes equations and is a well developed subject in modelling heat conducting fluids\cite{Z89,Z03,GLCS80,CD80}. The aim of this paper is to generalize the setting in \cite{OK16,OK17}, where $g$-B\'enard equations was considered, to time-fractional $g$-B\'enard equations.

This paper is organized as follows: In the first section, we will recall some concept and notations related to fractional calculus. Section 2, is devoted to the problem statement. Section 3 will be dedicated to the proof of the existence and uniqueness of weak solutions to time-fractional $g$-Navier-Stokes equations.

\section{Preliminaries on fractional calculus}

In this section, we provide some notations and preliminary results concerning fractional calculus. For this purpose, assume $X$ to be a Banach space. Let $\alpha\in (0,1]$ and let $k_\alpha$ denote the Riemann-Liouville kernel
\[
k_\alpha(t)=\frac{t^{\alpha-1}}{\Gamma(\alpha)}
\]

For a function $v:[0,T]\rightarrow X$, we give the following definitions of derivatives and integrals:
\begin{itemize}
\item[(1)] The left Riemann-Liouville integral of $v$ is defined by
\[
I_t^\alpha v(t)=\int_0^t k_\alpha(t-s)v(s)\,ds,\quad t>0
\]provided the integral is point-wise defined on $[0,+\infty[$
\item[(2)] The right Riemann-Liouville integral of $v$ is defined by
\[
\, I_{t,T}^\alpha v(t)=\int_t^T k_\alpha(t-s)v(s)\,ds,\quad t>0
\]provided the integral is point-wise defined on $[0,+\infty[$
\item[(3)]  The left Caputo fractional derivative of order $\alpha$ of $v$, is defined by
\[
\,D_t^\alpha v(t)= \int_0^t k_{1-\alpha}(t-s)\frac{d}{ds}v(s)\,ds
\]
\item[(4)] The right Riemann-Liouville fractional derivative of order $\alpha$ of $v$ is defined by
\[
\,D_{t,T}^\alpha v(t)=-\frac{d}{dt} \int_t^T k_{1-\alpha}(t-s)v(s)\,ds
\]
\item[(5)] The Liouville-Weyl fractional integral on the real axis for functions $v:\mathbb R\rightarrow X$ is defined as follows
\[
\,I_{-,t}^\alpha v(t)=\int_{-\infty}^t k_{\alpha}(t-s)v(s)\,ds
\]
\item[(6)] The Caputo fractional derivative on the real axis for functions $v:\mathbb R\rightarrow X$ is defined as follows
\[
\,D_{-,t}^\alpha v(t)= \,I_{-,t}^{1-\alpha}\frac{d}{dt} v(t)
\]
\end{itemize}

Note that the notation $\partial_t^\alpha$ stands for Caputo fractional partial derivative, i.e. when functions have another argument than time. We have the following fractional integration by parts formula; see, e.g. \cite{A07}
\begin{align}\label{partial}
\int_0^T(\partial_t^\alpha u(t),\psi(t))dt&=\int_0^T(u(t),D_{t,T}^\alpha\psi(t))dt+(u(t),I_{t,T}^{1-\alpha}\psi(t))|_0^T\\\nonumber
&=\int_0^T(u(t),D_{t,T}^\alpha\psi(t))dt-(u(0),I_T^{1-\alpha}\psi(t))
\end{align}since for $\psi\in C_0^\infty([0,T],X)$ one have $\displaystyle\lim_{t\to T}\,I_{t,T}^{1-\alpha}\psi(t)=0$.

To pass from weak convergence to strong convergence we will need a compactness result. Let $X_0$, $X$, $X_1$ be Hilbert spaces with $ X_0 \xhookrightarrow{} X \xhookrightarrow{} X_1$ being continuous and $X_0\xhookrightarrow{} X$ being compact. Assume that $v:\mathbb R\rightarrow X_1$ and denote by $\widehat v$ its Fourier transform
\[
\widehat v(\tau)=\int_{-\infty}^{+\infty}e^{-2i\pi t \tau}v(t)\,dt
\]
We have for $\gamma>0$
\[
\widehat{D_t^\gamma v}(\tau)=(2i\pi\tau)^\gamma \widehat v(\tau)
\]
For a given $0<\gamma<1$, we introduce the following space
\[
W^\gamma(\mathbb R,X_0,X_1)=\left\{v\in L^2(\mathbb R,X_0):\, D_t^\gamma v\in L^2(\mathbb R,X_1)\right\}
\]
Clearly, it is a Hilbert space for the norm
\[
\|v\|_\gamma=\left(\|v\|^2_{L^2(\mathbb R,X_0)}+\|\,|\tau|^\gamma\widehat v\|^2_{L^2(\mathbb R,X_1)}\right)^{1/2}
\]
 For any set $K\subset \mathbb R$, we associate with it the subspace $W_K^\gamma\subset W^\gamma$ defined as
 \[
 W_K^\gamma(\mathbb R,X_0,X_1)=\{v\in W^\gamma(\mathbb R,X_0,X_1):\,\mathrm{support }\, u\subset K\}
 \]

By similar discussion as in the proof of Theorem 2.2 in Temam \cite{T84}(see also Theorem 2.1 in \cite{ZP17}), it is clear that $W_K^\gamma(\mathbb R,X_0,X_1)\xhookrightarrow{} L^2(\mathbb R,X)$ is compact for any bounded set $K$ and any $\gamma>0$.

As a particular situation of the compactness result discussed above, let $H$, $V$ be two Hilbert spaces endowed with the scalar product $(.,.)_H$ and $(.,.)_V$ and the norms $|.|_H$ and $\|.\|_V$, respectively. Denote by $\langle.,.\rangle$ the dual pairing between $V$ and $V'$, the dual of $V$. Moreover assume that $V \xhookrightarrow{} H  \xhookrightarrow{} V'$ continuously and compactly and note that the space
\[
W^\gamma(0,T;V,V')=\left\{v\in L^2(0,T;V):\, \partial_t^\gamma v\in L^2(0,T;V')\right\}
\]is compactly embedded in $L^2(0,T;H)$. Similarly to Lemma 2.1 in \cite{ZP17}, we have
\[
\partial_t^\gamma(u(t),v)_V=\langle\partial_t^\gamma u(t),v\rangle
\] for $u\in W^\gamma(0,T;V,V')$ and $v\in H$. Moreover, for a derivable function $v:[0,T]\rightarrow V$ we have from \cite{A10} that

\[
(v(t),D_t^\gamma v(t))_{H}\geq\frac{1}{2}\,D_t^\gamma|v(t)|^2
\]

We end this section by the following important result
\begin{lemma}

Suppose that a nonnegative function satisfies
\[
_{~0}^{~C}D_t^\gamma v(t)+ c_1 v(t)\leq c_2(t)
\]for $c_1>0$ and $c_2$ a nonnegative integrable function for $t\in[0,T]$. Then
\[
v(t)\leq v(0)+\frac{1}{\Gamma(\gamma)}\int_0^t(t-s)^{\gamma-1} c_2(s)\,ds
\]

\end{lemma}

For more details about fractional calculus we refer to the monographs \cite{KST06,Z14,Z16}.

\section{ Problem Statement}

We introduce the usual notation used in the context of the mathematical theory of Navier-Stokes equations \cite{T84}. Let $\Omega_g = \Omega_2 \times (0, g) = (0, 1)\times (0, 1)\times (0, g)$ where
$g = g(y_1, y_2)$ is a smooth function defined on $\Omega_2$. In addition we assume that

\begin{equation}\label{2.1}
  \begin{array}{c }
\displaystyle  0<m_0< g(y_1,y_2)\leq M_0,~~~~~~\text{for all}~~ (y_1, y_2)\in \Omega_2
  \\ [0.2cm]
  \displaystyle |\nabla g|_{\infty}=\displaystyle\sup_{\Omega_2}| \nabla g  |<\infty,~~~~~~~~g\in C^{\infty}_{per}(\Omega_2).
\end{array}
\end{equation}

Let $L^2(\Omega, g)$ denotes the Hilbert space, of weighted Sobolev spaces type, with the inner product \[\langle u,v\rangle_g = \int_\Omega(u \cdot v)gdx\] and the induced norm $|u|^2_g = \langle u, u\rangle_g$. Similarly, we can define the weithed Sobolev space $H^1(\Omega, g)$ equipped with the norm

\[ |u|^2_{H^1(\Omega,g)}=\langle u,u \rangle_g + \displaystyle \sum^n_{i=1} \langle \frac{\partial u}{\partial x_i} ,\frac{\partial u}{\partial x_i} \rangle_g.  \]
Moreover, we will need the following spaces:

\begin{align*}
\mathcal{V}_1 & = \, \{ u\in (C^{\infty}_{per}(\Omega))^n~~: ~~ \nabla\cdot(gu)=0,~~~\displaystyle\int_{\Omega}udx=0~~~\text{on}~~~\Omega      \}\\
H_g & =  \text{ the closure of } \mathcal{V}_1 \text{ on}~~~ L^2(\Omega,g)\\
 V_g & =   \text{ the closure of }\mathcal{V}_1 \text{ on }~~~ H^1(\Omega,g) \\
 V'_g & =    \text{  the dual space of  }~~~V_g \\
  \mathcal{V}_2 & = \, \{ \varphi\in C^{\infty}_{per}(\Omega)~~:~~ \displaystyle\int_{\Omega}\varphi dx=0    \}\\
 W_g & =  \text{ the closure of } \mathcal{V}_2 \text{ on}~~~ H^1(\Omega,g) \\
  W'_g & =  \text{  the dual space of  }~~~W_g\\
    Q & =  \text{  the closure of}~~~ \{\nabla\varphi~~:~~\varphi\in C^1_{per}(\overline{\Omega},R)  \} ~~\text{in}~~L^2(\Omega)
\end{align*}
where $H_g$ is endowed with the inner product and the norm in $L^2(\Omega g)$. In addition, the spaces  $V_g$ and $Wg$ are endowed with the inner product and the norm in $H^1(\Omega, g)$. Let us also remark that the inclusions
\[V_g\subset  H_g = H'_g\subset V'_g,\] \[W_g \subset L^2(\Omega, g) \subset W'_g,\]
are dense and continuous \cite{M92,R01}. By the Riesz representation theorem, it is possible to write
\[ \langle f,u \rangle_g=(f,u)_g,~~~~~~~~ \forall f\in H_g,~\forall u\in V_g.      \]

Let us now define the orthogonal projection $P_g$ as $P_g : L^2_{ per}(\Omega, g)\rightarrow  H_g$. It is clear that $Q \subseteq H^{\perp}_{g}$ . Similarly, we define $\widetilde{P}_g$ as $\widetilde{P}_g :L^2_{ per}(\Omega, g)\rightarrow  W_g$. By taking into account the following equality\cite{R01}:
\[  \displaystyle -\frac{1}{g} ( \nabla\cdot g\nabla u )  =-\Delta u - \frac{1}{g}( \nabla g\cdot \nabla  )u\]we define the $g-$Laplace operator and $g$-Stokes operator as follows
\[-\Delta_gu=-\displaystyle \frac{1}{g}( \nabla\cdot g\nabla u)  \] and
\[A_gu = P_g[-\Delta_gu]\] respectively. We have the following result \cite{R01}:

\begin{proposition} For the g-Stokes operator $A_g$, the followings hold:

\begin{enumerate}
  \item[(1)] The g-Stokes operator $A_g$ is a positive, self-adjoint operator with compact inverse, where the domain of $A_g$ is
  $D(A_g) = V_g \cap H^2(\Omega, g).$
  \item[(2)]  There exist countable eigenvalues of $A_g$ satisfying
  \[     0< \displaystyle \frac{4\pi^2m_0}{M_0}\leq  \lambda_1  \leq \lambda_2 \leq\lambda_3\leq\cdot\cdot\cdot   \]
where $\lambda_1$ is the smallest eigenvalue of $A_g$. In addition, there exist the corresponding collection of eigenfunctions $\{w_i\}_{i\in\mathbb{N}} $ forms an orthonormal basis for
$H_g$.
\end{enumerate}
\end{proposition}

The operators $A_g$ and $P_g$ are clearly self-adjoint, then by using integration by parts we have:
  \begin{align*}
  \langle A_g u,u\rangle_g= &\langle P_g[ \displaystyle -\frac{1}{g}(\nabla\cdot g \nabla )u],u\rangle_g\\
    =&\displaystyle \int_{\Omega}  (\nabla u\cdot\nabla u )_gdx \\
    = &  \langle \nabla u\cdot\nabla u \rangle_g
 \end{align*}
It then follows that for $u \in V_g$ we can write $  | A^{1/2}u|_g=  | \nabla u |_g=\| u \|_g   $. On the other hand,, since the functional
\[  \tau\in W_g\rightarrow ( \nabla \theta,\nabla\tau)_g\in \mathbb{R}    \]
is a continuous linear mapping on $W_g$, we can define a continuous linear mapping $\widetilde{A}_g$ on $W'_g$ such that
\[    \forall \tau \in W_g,~~~~~\langle \widetilde{A}_g,\tau \rangle_g = ( \nabla \theta  , \nabla \tau )_g \]
for all $\theta\in  W_g$. For $u$, $v$ and $w$ laying in an appropriate subspaces of $L^2_{per}(\Omega, g)$, we can define the bilinear operator \[B_g(u, v) = P_g[(u \cdot \nabla)v]\] and the trilinear form
\[  b_g(u,v,w)=\displaystyle\sum^{n}_{i,j=1}\int_{\Omega}   u_i\frac{\partial v_j}{\partial x_i} w_j gdx \]
As a consequence, one obtains $b_g(u, v, w) =
-b_g(u, w, v)$ which implies that $b_g(u, v,v) = 0$. MOreover, we have the following inequality satisfied by $b_g$(see \cite{T84,W09}),
\begin{equation}\label{2.7}
| b_g(u,v,w)|_g\leq   c| u|^{1/2}_{g}  \| u\|^{1/2}_{g}  | v|_{g}   | w|^{1/2}_{g}  \| w\|^{1/2}_{g},~~~~~~\forall u,v,w\in V_g
\end{equation}

Similarly, for $u \in V_g$ and $\theta, \tau\in W_g$ we define $\widetilde{B}_g(u,\theta) = \widetilde{P}_g[(u \cdot \nabla)\theta]$ and
\[  \widetilde{b}_g(u,\theta,\tau)=\displaystyle\sum^{n}_{i,j=1}\int_{\Omega}   u_i(x)\frac{\partial \theta(x)}{\partial x_j} \tau(x) gdx \]

We denote the operators
$  C_gu = P_g[\frac{1}{g} (\nabla g \cdot \nabla)u]$  and $\widetilde{C}_g \theta = \widetilde{P}_g[\frac{1}{g} (\nabla g \cdot \nabla)\theta]$ such that
\begin{align*}
 \langle  C_gu,v \rangle_g  =&b_g (\displaystyle \frac{\nabla g}{g},u,v),\\
 \langle \widetilde{C}_g \theta,\tau  \rangle_g=&\widetilde{b}_g (\displaystyle \frac{\nabla g}{g},\theta,\tau  )
\end{align*}
Finally, let $\widetilde{D}_g\theta= \widetilde{P}_g[\displaystyle \frac{\nabla g}{g}\theta ]$ such that
\[ \langle  \widetilde{D}_g\theta ,\tau  \rangle_g=-\widetilde{b}_g (\displaystyle \frac{\nabla g}{g},\theta,\tau  ) - \widetilde{b}_g (\displaystyle \frac{\nabla g}{g},\theta,\tau  ). \]
We can now rewrite the system of $g$-B\'enard equations in the following abstract time-fractional evolutionary equations
\begin{equation}\label{benard}
\begin{array}{ccc}
&\partial^{\alpha}_{t}u+ B_g(u,u)+\nu A_gu+\nu C_gu= \xi\theta+f_1\\[0.2cm]
&\partial^{\alpha}_{t} \theta+ \widetilde{B}_g(u,\theta)+\kappa \widetilde{A}_g\theta-\kappa \widetilde{C}_g\theta-   \kappa\widetilde{D}_g\theta= f_2\\[0.2cm]
&u(x,0)=u_0(x),~~~~~~~~\theta(x,0)=\theta_0(x)
\end{array}
\end{equation}

We give the following two lemmas and for the proofs of these lemmas we refer to \cite{BR04,L69}.

\begin{lemma}\label{lem1} For $n = 2$ there exists a positive constant $c$ such that

\[ |u|_{L^4(\Omega,g)} \leq c|u|^{1/2}_{g}  |\nabla u|^{1/2}_{g} ,~~~~~~~~ \forall u\in H^1( \Omega , g ) . \]
\end{lemma}

\begin{lemma}\label{lem2} For  $   u\in L^2(  0,T,V_g)$ we have:
\[B_g(u,u)(t)\in L^1(0,T,V'_g)  \text{  and  } C_gu(t)\in L^2(0,T,H_g).      \]
\end{lemma}

\section{    Existence of weak solutions   }

\begin{definition} A pair of functions $ \{ u,\theta\}$ is called a weak solution of the system \eqref{benard} if $u\in  L^2(0, T; V_g)$ and $\theta \in L^2(0, T; W_g)$ satisfy the following equations

\begin{equation}\label{3.1}
\begin{array}{ccc}
&\partial^{\alpha}_{t}(u,v)_g+ b_g(u,u,v)+\nu (\nabla u,\nabla v  )_g+\nu (C_gu,v)_g= (\xi\theta,v)_g+(f_1,v)_g\\[0.2cm]
&\partial^{\alpha}_{t} (\theta,\tau)_g+ \widetilde{b}_g(u,\theta,\tau)+\kappa(  \nabla\theta,\nabla\tau )_g+\kappa \widetilde{b}_g
( \displaystyle \frac{\nabla g}{g},\tau,\theta)=( f_2,v)_g
\end{array}
\end{equation}

for all $v_2 \in V_g$ and $\tau\in W_g.$
\end{definition}
\begin{theorem}If $f_1\in  L^{\frac{2}{\alpha_1}}(0, T; L^2(\Omega, g))$ and $f_2\in  L^{\frac{2}{\alpha_2}}(0, T; L^2(\Omega, g))$ $(\alpha_1,\,\alpha_2<\alpha)$, $u_0 \in H_g$, $\theta_0 \in L^2(\Omega, g)$ and $g$ is a smooth function satisfying the conditions given in (\ref{2.1}) defined on $\Omega_2$ then, there exist a unique weak solution $\{ u,\theta \}  $ of the system \eqref{benard} satisfying the periodic boundary conditions.
\end{theorem}

For the proof of the theorem, we shall use the standard Feado- Galerkin method.

\begin{proof}

Since $V_g$ is separable and $\mathcal{V}_1$ is dense in $V_g$, there exists a sequence $\{u_i\}_{i\in \mathbb{N}}$ which forms a
complete orthonormal system in $H_g$ and a base for $V_g$. Similarly, there exists a sequence
$\{\theta_i\}_{i\in \mathbb{N}}$ which forms a complete orthonormal system in $L^2(\Omega, g)$ and a base for $W_g$. Let $m$ be an arbitrary but fixed positive integer. For each $m$, we define an approximate
solution $\{u^{(m)}(t), \theta^{(m)}(t)\}$ of \eqref{benard} for $1 \leq k \leq m$ and $t \in [0 ,T]$ in the form,

\begin{equation}\label{3.3}
\begin{array}{lr}
 \displaystyle u^{(m)}(t)=\sum^m_{j=1} f_j^{(m)}(t)u_j,~~~~~~~& \displaystyle\theta^{(m)}(t)=\sum^m_{j=1} g_j^{(m)}(t)\theta_j
\end{array}
\end{equation}
and we consider the following approximate problem

\begin{equation}\label{3.4}
\begin{array}{l}
\partial^{\alpha}_{t}(u^{(m)},u_k)_g+ b_g(u^{(m)},u^{(m)},u_k)+\nu ((u^{(m)},u_k ) )_g
+\nu b_g( \displaystyle \frac{\nabla g}{g}   ,u^{(m)},u_k)\\
\qquad= (\xi\theta^{(m)},u_k    )_g+(f_1,u_k)_g,
\end{array}
\end{equation}

\begin{equation}\label{3.5}
\begin{array}{l}
\partial^{\alpha}_{t} (\theta^{(m)} ,\theta_k)_g+ \widetilde{b}_g(u^{(m)},\theta^{(m)},\theta_k)+\kappa(( \theta^{(m)},\theta_k ))_g
+\kappa \widetilde{b}_g( \displaystyle \frac{\nabla g}{g},\theta_k,\theta^{(m)})\\
\qquad=( f_2,\theta_k)_g,
\end{array}
\end{equation}and
\begin{equation}\label{3.6}
\begin{array}{lr}
 \displaystyle u^{(m)}(0)=u_{m_0}=\sum^m_{j=1} (a_0,u_j)u_j,~~~~~~~& \displaystyle\theta^{(m)}(0)=\theta_{m_0}=\sum^m_{j=1} (\tau_0,\theta_j)\theta_j
\end{array}
\end{equation}
This system forms a nonlinear fractional order system of ordinary differential equations for the functions $f_j^{(m)}(t)$ and $g_j^{(m)}(t)$ and has a maximal solution on some interval $[0, T ]$ ( cf. \cite{BA20}). We multiply (\ref{3.4}) and (\ref{3.5}) by $f_j^{(m)}(t)$ and $g_j^{(m)}(t)$ respectively, and add these equations for $k = 1,\cdot\cdot\cdot, m.$ Taking into account $b_g(u^{(m)}, u^{(m)}, u^{(m)}) = 0$ and $\widetilde{b}g(u^{(m)}, \theta^{(m)},\theta^{(m)}) =0$ we get;

\begin{equation}\label{3.7}
\begin{array}{r}
(D^{\alpha}_{t}u^{(m)},u^{(m)})_g+ \nu\| u^{(m)} (t) \|^2_g +\nu b_g( \displaystyle \frac{\nabla g}{g}   ,u^{(m)}(t),u^{(m)}(t))
\\ = (\xi\theta^{(m)},u^{(m)}(t)    )_g+(f_1,u^{(m)}(t))
\end{array}
\end{equation}
and
\begin{equation}\label{3.8}
\begin{array}{r}
 (D^{\alpha}_{t}\theta^{(m)}(t) ,\theta^{(m)}(t))_g+  \kappa\| \theta^{(m)}(t)  \|^2_g   +\kappa \widetilde{b}_g( \displaystyle \frac{\nabla g}{g},\theta^{(m)}(t),\theta^{(m)}(t))=( f_2,\theta^{(m)}(t))_g
\end{array}
\end{equation}

Using Schwarz and Young inequalities in (\ref{3.7}) and (\ref{3.8})

\begin{align*}
&D^{\alpha}_t  |  u^{(m)}(t) |^2_{g} +\nu \| u^{(m)}(t) \|^2_{g}  \leq   \displaystyle \frac{ M_0 |\xi|^2_{\infty}}{\pi^2 m_0\nu} |  \theta^{(m)}(t)   |^2_g                +\frac{4}{\nu}\| f_1(t) \|^2_{V'_g}+\frac{   2\nu |\nabla g|^2_{\infty}}{m^2_0}|   u^{(m)}(t) |^2_g\\
&D^{\alpha}_t  |  \theta^{(m)}(t) |^2_{g} +\kappa \| \theta^{(m)}(t) \|^2_{g}  \leq    \displaystyle \frac{ 2}{\kappa}\| f_2(t) \|^2_{W'_g}
+\frac{2\kappa|\nabla g|^2_{\infty}}{m^2_0}|\theta^{(m)}(t)|^2_g
\end{align*}
By using the fact that  $\displaystyle |\nabla g |^2_{\infty}< \frac{\pi^2 m^3_0}{M_0}$ and noting $\nu'=\nu (1- \frac{ M_0|\nabla g|^2_{\infty}}{2\pi^2 m^3_0} )$, $\kappa'=\kappa(1- \frac{ M_0|\nabla g|^2_{\infty}}{2\pi^2 m^3_0}      )$ and $c'= \frac{ M^2_0\|\xi\|^2_{\infty}}{4\pi^4 m^2_0}$, we get the inequalities
\begin{equation}\label{3.9}
\begin{array}{l}
D^{\alpha}_t  |  u^{(m)}(t) |^2_{g} +\nu' \| u^{(m)}(t) \|^2_{g}  \leq   \displaystyle \frac{ c'}{\nu} \|  \theta^{m}(t)   \|^2_g                +\frac{4}{\nu}\| f_1(t) \|^2_{V'_g}
\end{array}
\end{equation}
and
\begin{equation}\label{3.10}
\begin{array}{l}
 D^{\alpha}_t  |  \theta^{(m)}(t) |^2_{g} +\kappa' \| \theta^{(m)}(t) \|^2_{g}  \leq    \displaystyle \frac{ 2}{\kappa}\| f_2(t) \|^2_{W'_g}  \end{array}
\end{equation}
Integrating  \eqref{3.10} from $0$ to $T$, in the fractional sense, we obtain

\begin{align*}
|\theta^{(m)}(t)|_g^2 +&\frac{\kappa'}{\Gamma(\alpha)}\int_0^t (t-s)^{\alpha-1}\|\theta^{(m)}(s)\|_g^2\, ds\\\leq&\, |\theta_{0m}|_g^2+\frac{2}{\kappa\Gamma(\alpha)}\int_0^t(t-s)^{\alpha-1}\| f_2(s)\|_{W_g'}^2\,ds\\
\leq &\, |\theta_{0m}|_g^2+ \frac{2}{\kappa\Gamma(\alpha)}\int_0^t\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+\frac{2}{\kappa\Gamma(\alpha)}\int_0^t(t-s)^{\frac{\alpha-1}{1-\alpha_2}}\,ds\\
\leq &\, |\theta_{0m}|_g^2+ \frac{2}{\kappa\Gamma(\alpha)}\int_0^T\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+C_2
\end{align*}
where $b_2=\frac{\alpha-1}{1-\alpha_2}$ and $C_2=\frac{2T^{1+b_2}}{\kappa(1+b_2)\Gamma(\alpha)}$. It follows that
\begin{equation}
\int_0^t (t-s)^{\alpha-1}\|\theta^{(m)}(s)\|_g^2\, ds\leq\, \frac{\Gamma(\alpha)}{\kappa'}|\theta_{0m}|_g^2+\frac{2}{\kappa\kappa'}\int_0^T\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+\frac{\Gamma(\alpha)}{\kappa'}C_2
\end{equation}
On the other hand, integrating  \eqref{3.9} from $0$ to $T$, in the fractional sense, we obtain
\begin{align*}
|u^{(m)}(t)|_g^2 &+\frac{\nu'}{\Gamma(\alpha)}\int_0^t (t-s)^{\alpha-1}\|u^{(m)}(s)\|_g^2\, ds\\
\leq&\, |u_{0m}|_g^2+\frac{c'}{\nu\Gamma(\alpha)}\int_0^t (t-s)^{\alpha-1}\|\theta^{(m)}(s)\|_g^2\, ds+\frac{4}{\nu\Gamma(\alpha)}\int_0^t(t-s)^{\alpha-1}\| f_1(s)\|_{V_g'}^2\,ds\\
\leq &\,|u_{0m}|_g^2+\frac{c'}{\nu\kappa'}|\theta_{0m}|_g^2+ \frac{2c'}{\nu\kappa\kappa'\Gamma(\alpha)}\int_0^t\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+\frac{c'}{\nu\kappa'}C_2\\&+\frac{4}{\nu\Gamma(\alpha)}\int_0^t\|f_1(s)\|_{V_g'}^{2/\alpha_1}\,ds+\frac{4}{\nu\Gamma(\alpha)}\int_0^t(t-s)^{\frac{\alpha-1}{1-\alpha_1}}\,ds\\
\leq &\,|u_{0m}|_g^2+\frac{c'}{\nu\kappa'}|\theta_{0m}|_g^2+ \frac{2c'}{\nu\kappa\kappa'\Gamma(\alpha)}\int_0^t\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+\frac{4}{\nu\Gamma(\alpha)}\int_0^t\|f_1(s)\|_{V_g'}^{2/\alpha_1}\,ds\\
&+C_1
\end{align*}
where $b_1=\frac{\alpha-1}{1-\alpha_1}$ and $C_1=\frac{c'}{\nu\kappa'}C_2+\frac{4T^{1+b_1}}{\nu(1+b_1)\Gamma(\alpha)}$. By using the fact that
\begin{equation}
\int_0^t (t-s)^{\alpha-1}\|u^{(m)}(s)\|_g^2\, ds\geq T^{\alpha-1}\int_0^t\|u^{(m)}(s)\|_g^2\, ds
\end{equation}and similarly
\begin{equation}
\int_0^t (t-s)^{\alpha-1}\|\theta^{(m)}(s)\|_g^2\, ds\geq T^{\alpha-1}\int_0^t\|\theta^{(m)}(s)\|_g^2\, ds
\end{equation}it follows that
\begin{align}
|u^{(m)}(t)|_g^2 +\frac{\nu' T^{\alpha-1}}{\Gamma(\alpha)}&\int_0^t\|u^{(m)}(s)\|_g^2\, ds\leq
|u_{0m}|_g^2+\frac{c'}{\nu\kappa'}|\theta_{0m}|_g^2\\
&+ \frac{2c'}{\nu\kappa\kappa'\Gamma(\alpha)}\int_0^T\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds+\frac{4}{\nu\Gamma(\alpha)}\int_0^T\|f_1(s)\|_{V_g'}^{2/\alpha_1}\,ds +C_1\nonumber
\end{align}

\begin{equation}
|\theta^{(m)}(t)|_g^2 +\frac{\kappa'T^{\alpha-1}}{\Gamma(\alpha)}\int_0^t\|\theta^{(m)}(s)\|_g^2\, ds\leq |\theta_{0m}|_g^2+ \frac{2}{\kappa\Gamma(\alpha)}\int_0^T\|f_2(s)\|_{V_g'}^{2/\alpha_2}\,ds+C_2
\end{equation}
Consequently
\begin{align}
\displaystyle\sup_{t\in[0,T]}|u^{(m)}(t)|_g^2 \leq
|u_{0m}|_g^2+&\frac{c'}{\nu\kappa'}|\theta_{0m}|_g^2+ \frac{2c'}{\nu\kappa\kappa'\Gamma(\alpha)}\int_0^T\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds\\&+\frac{4}{\nu\Gamma(\alpha)}\int_0^T\|f_1(s)\|_{V_g'}^{2/\alpha_1}\,ds +C_1\nonumber
\end{align}

\begin{equation}\label{3.12}
 \displaystyle \sup_{t\in[0,T]}| \theta^{(m)}(t) |^2_g
\leq |\theta_{0m}|_g^2+ \frac{2}{\kappa\Gamma(\alpha)}\int_0^T\|f_2(s)\|_{V_g'}^{2/\alpha_2}\,ds+C_2
\end{equation}
which imply that the sequences $\{u^{(m)}\}_m$ and $\{\theta^{(m)}\}_m$ remain in a bounded set of
$L^{\infty}(0, T ; H_g)$ and $L^{\infty}(0, T ; L^2(\Omega, g))$ respectively. Moreover, for $t=T$, one obtain
\begin{align}\label{u}
\int_0^T\|u^{(m)}(s)\|_g^2\, ds\leq&
\frac{\Gamma(\alpha)}{\nu' T^{\alpha-1}}|u_{0m}|_g^2+\frac{c'}{\nu\nu'\kappa'T^{\alpha-1}}|\theta_{0m}|_g^2+ \frac{2c'}{\nu\nu'\kappa\kappa'T^{\alpha-1}}\int_0^T\|f_2(s)\|_{W_g'}^{2/\alpha_2}\,ds\\
&+\frac{4}{\nu\nu'T^{\alpha-1}}\int_0^T\|f_1(s)\|_{V_g'}^{2/\alpha_1}\,ds +\frac{\Gamma(\alpha)}{\nu' T^{\alpha-1}}C_1\nonumber
\end{align}

\begin{equation}\label{theta}
\int_0^T\|\theta^{(m)}(s)\|_g^2\, ds\leq \frac{\Gamma(\alpha)}{\kappa'T^{\alpha-1}}|\theta_{0m}|_g^2+ \frac{2}{\kappa\kappa'T^{\alpha-1}}\int_0^T\|f_2(s)\|_{V_g'}^{2/\alpha_2}\,ds+\frac{\Gamma(\alpha)}{\kappa'T^{\alpha-1}}C_2
\end{equation}
which imply that the sequences  $\{u^{(m)}\}_m$ and  $\{\theta^{(m)}\}_m$ remain in a bounded set of
$L^2(0, T ; V_g)$ and $L^2(0, T ; W_g)$ respectively. Consequently, we can assert the existence of elements
$u \in  L^2(0, T ; V_g) \cap L^{\infty}(0, T ; H_g)$ and $\theta \in L^2(0, T ; W_g) \cap L^{\infty}(0, T ; L^2(\Omega, g))$ and the subsequences  $\{u^{(m)}\}_m$ and  $\{\theta^{(m)}\}_m$ such that $u^{(m)} \rightarrow u \in
L^2(0, T ; V_g)$ and $\theta^{(m)} \rightarrow\theta \in L^2(0, T ; W_g)$ weakly and $u^{(m)}\rightarrow u \in L^{\infty}(0, T ; H_g)$ and
$\theta^{(m)} \rightarrow \theta\in  L^{\infty}(0, T ; L^2(\Omega, g))$ weak-star
 as $m \rightarrow\infty$.

Let $\widetilde{u}^{(m)} : \mathbb{R} \rightarrow V_g$
 and $\widetilde{\theta}^{(m)} : \mathbb{R} \rightarrow W_g$ defined as

\[   \widetilde{u}^{(m)}(t)=\left\{
  \begin{array}{rcl}
u^{(m)}(t),&0\leq t\leq T\\
0,& \text{otherwise}
    \end{array} \right.~~~~~\text{and}~~~~~~
      \widetilde{\theta}^{(m)}(t)=\left\{     \begin{array}{rcl}
\theta^{(m)}(t),&0\leq t\leq T\\
0,& \text{otherwise}
    \end{array} \right.     \]
and their Fourier transforms denoted by $\widehat{u}^{(m)}$ and  $\widehat{\theta}^{(m)}$, respectively. We show that the sequence $\{\tilde u^{(m)}\}_m$ remains bounded in $W^\gamma(\mathbb R, V_g, H_g)$ and the sequence $\{\tilde \theta^{(m)}\}_m$ remains bounded in $W^\gamma(\mathbb R, W_g, L^2(\Omega,g))$. To do so, we need to verify that
\begin{equation}\label{eqgamma1}
\int_{-\infty}^{+\infty}|\tau|^{2\gamma}|\widehat u^{(m)}(\tau)|^2\,d\tau\leq \mathrm{const.}\quad\text{for some }\gamma>0
\end{equation}and
\begin{equation}\label{eqgamma2}
\int_{-\infty}^{+\infty}|\tau|^{2\gamma}|\widehat \theta^{(m)}(\tau)|^2\,d\tau\leq \mathrm{const.}\quad\text{for some }\gamma>0
\end{equation}
In order to prove \eqref{eqgamma1} and \eqref{eqgamma2}, we observe that
\begin{equation}\label{Fmu}
(D_t^\alpha\tilde u^{(m)},u_k)_g=(\widetilde F^u_m,u_k)_g+(u_{m0},u_k)_g\,I_{-,t}^{1-\alpha}\delta_0-( u^{(m)}(T),u_k)_g \,I_{-,t}^{1-\alpha}\delta_T
\end{equation}
\begin{equation}\label{Fmtheta}
(D_t^\alpha\tilde \theta^{(m)},\theta_k)_g=(\widetilde F^\theta_m,\theta_k)_g+(\theta_{m0},\theta_k)_g\,I_{-,t}^{1-\alpha}\delta_0-( \theta^{(m)}(T),\theta_k)_g \,I_{-,t}^{1-\alpha}\delta_T
\end{equation}
where $\delta_0,\,\delta_T$ are Dirac distributions at $0$ and $T$ and $F_m^u$ and $F_m^\theta$ are defined by
\[
\begin{array}{rcl}
F_m^u=&   \xi \theta^{(m)}+f_1-B_g(u^{(m)},u^{(m)})-\nu A_g u^{(m)}-\nu C_gu^{(m)}\\
F_m^\theta=& f_2-\widetilde{B}_g(u^{(m)},\theta^{(m)}   )-\kappa \widetilde{A}_g\theta^{(m)}+\kappa\widetilde{C}_g\theta^{(m)}+\kappa \widetilde{D}_g\theta^{(m)}
\end{array}
\]for $k = 1,\cdot,\cdot\cdot, m$. Here $\widetilde F_m$ is defined as usual by

\begin{equation}
 \widetilde F_m(t)=
 \left \{
   \begin{array}{r c l}
    F_m(t),&\quad 0\leq t\leq T   \\
    0,&\quad \text{ otherwise}
   \end{array}
   \right .
\end{equation}

Indeed, it is classical that since $\widetilde u^{(m)}$ and $\widetilde \theta^{(m)}$ have two discontinuities at $0$ and $T$, the Caputo derivative of $\widetilde u^{(m)}$ is given by

\begin{align}
D_{-,t}^\alpha\tilde u^{(m)}&= \,I_{-,t}^{1-\alpha}\left(\frac{d}{dt}\tilde u^{(m)}\right)\\
&=\,I_{-,t}^{1-\alpha}\left(\frac{d}{dt} u^{(m)}+u^{(m)}(0)\delta_0-u^{(m)}(T)\delta_T\right)\\
&=\,D_t^\alpha u^{(m)}+I_{-,t}^{1-\alpha}\left(u^{(m)}(0)\delta_0-u^{(m)}(T)\delta_T\right)
\end{align}and the one of $\widetilde \theta^{(m)}$ is given by
\begin{equation}
D_{-,t}^\alpha\widetilde \theta^{(m)}=\,D_t^\alpha \theta^{(m)}+I_{-,t}^{1-\alpha}\left(\theta^{(m)}(0)\delta_0-\theta^{(m)}(T)\delta_T\right)
\end{equation}

By the Fourier transform, \eqref{Fmu} and \eqref{Fmtheta} yield to
\begin{align}\label{fourieru}
(2i\pi\tau)^\alpha(\widehat u^{(m)},u_k)_g=&(\widehat{F}^u_m,u_k)_g+(u_{m0},u_k)_g(2i\pi\tau)^{\alpha-1}\\
&\qquad\qquad-(u^{(m)}(T),u_k)_g(2i\pi\tau)^{\alpha-1}e^{-2i\pi T\tau}
\end{align}
\begin{align}\label{fouriertheta}
(2i\pi\tau)^\alpha(\widehat \theta^{(m)},\theta_k)_g=&(\widehat{F}^\theta_m,\theta_k)_g+(\theta_{m0},\theta_k)_g(2i\pi\tau)^{\alpha-1}\\
&\qquad\qquad-(\theta^{(m)}(T),\theta_k)_g(2i\pi\tau)^{\alpha-1}e^{-2i\pi T\tau}
\end{align}here $\widehat u^{(m)}$ and $\widehat{F}_m$ denote the Fourier transforms of $\tilde u^{(m)}$ and $\widetilde{F}_m$, respectively. We multiply \eqref{fourieru} and \eqref{fouriertheta} by $\widehat{f}_j^{(m)}$ and $\widehat{g}_j^{(m)}$ respectively and add
these equations for $k= 1,\dots, m$ to get
\begin{align}
(2i\pi\tau)^\alpha|\widehat u^{(m)}(\tau)|_g^2=&(\widehat{F}^u_m(\tau),\widehat u^{(m)}(\tau))_g+(u_{m0},\widehat u^{(m)}(\tau))_g(2i\pi\tau)^{\alpha-1}\\
&\qquad\qquad-(u^{(m)}(T),\widehat u^{(m)}(\tau))_g(2i\pi\tau)^{\alpha-1}e^{-2i\pi T\tau}
\end{align}
\begin{align}
(2i\pi\tau)^\alpha|\widehat \theta^{(m)}(\tau)|_g^2=&(\widehat{F}^\theta_m(\tau),\widehat \theta^{(m)}(\tau))_g+(\theta_{m0},\widehat \theta^{(m)}(\tau))_g(2i\pi\tau)^{\alpha-1}\\
&\qquad\qquad-(\theta^{(m)}(T),\widehat \theta^{(m)}(\tau))_g(2i\pi\tau)^{\alpha-1}e^{-2i\pi T\tau}
\end{align}
Since the integrals on the right hand side of the inequalities
\begin{align}
\displaystyle \int^T_0 \|F_m^u(t) \|_{V'_g}dt\leq \int^T_0 c( |\xi|_{\infty} \| \theta^{(m)}(t) \|_g&+\|f_1(t)\|_{V'_g}+| u^{(m)}(t)|_g \|u^{(m)} \|_g\\
&+ \|u^{(m)}(t)  \|_g+| \nabla g |_{\infty}\|u^{(m)}(t)  \|_g)dt\nonumber
\end{align}
\begin{align}
\displaystyle \int^T_0 \|F_m^\theta(t) \|_{W'_g}dt\leq \int^T_0 c'(  \| f_2(t) \|_{W'_g}+& |u^{(m)}(T)|_g \|\theta^{(m)} (t) \|_g+ \|\theta^{(m)} (t) \|_g\\
&  +|\nabla g|_{\infty} \| \theta^{(m)}(t) \|_g+|\Delta g|_{\infty}\| \theta^{(m)}(t) \|_g)dt\nonumber
\end{align}remains bounded, $ \|F_1(t) \|_{V'_g}$ and $\|F_2(t) \|_{W'_g}$ are bounded in $L^1(0, T ; V'_g)$ and $L^1(0, T ; W'_g)$ respectively. Therefore, for all $m$

\[
\begin{array}{ ccc}
\displaystyle \sup_{\tau\in\mathbb{R}} \| \widehat{F}^u_m(\tau) \|_{V'_g} \leq c_1&\text{and}&  \displaystyle \sup_{\tau\in\mathbb{R}} \| \widehat{F}^\theta_m(\tau) \|_{W'_g}\leq c_2.
\end{array}\]
Moreover, since $u^{(m)}(0)$,  $u^{(m)}(T),$  $\theta^{(m)}(0)$ and $\theta^{(m)}(T )$ are bounded, we get

\begin{align*}
    |\tau|^\alpha | \widetilde{u}^{(m)}(\tau)|^2_g      \leq &\, c_1\| u^{(m)} \|_{V_g}+c_2 |\tau|^{\alpha-1}\,|u^{(m)}  |_g\\
    \leq&\, c_3 \| u^{(m)} \|_{V_g}
    \end{align*}
\begin{align}
    |\tau|^\alpha | \widetilde{\theta}^{(m)}(\tau)|^2_g      \leq &\, c'_1\| \theta^{(m)} \|_{W_g}+c'_2 |\tau|^{\alpha-1}\,|\theta^{(m)}  |_g\\
    \leq &\,c_3 \| \theta^{(m)} \|_{W_g}
\end{align}
For $\gamma$ fixed, $\gamma<\alpha/4$, we observe that
\[
|\tau|^{2\gamma}\leq c(\gamma)\frac{1+|\tau|^\alpha}{1+|\tau|^{\alpha-2\gamma}}
\]
Then we can write
\begin{align*}
\int_{-\infty}^{+\infty}|\tau|^{2\gamma}|\widehat u^{(m)}(\tau)|^2_g\leq & c_5(\gamma) \int_{-\infty}^{+\infty}\frac{1+|\tau|^\alpha}{1+|\tau|^{\alpha-2\gamma}}|\widehat u^{(m)}(\tau)|^2_g\,d\tau\\
\leq & c_6(\gamma) \int_{-\infty}^{+\infty}\frac{1}{1+|\tau|^{\alpha-2\gamma}}\|\widehat u^{(m)}(\tau)\|^2_{V_g}\,d\tau\\
&\qquad\qquad+c_7(\gamma) \int_{-\infty}^{+\infty}\frac{|\tau|^{\alpha-1}}{1+|\tau|^{\alpha-2\gamma}}\|\widehat u^{(m)}(\tau)\|^2_{V_g}\,d\tau
\end{align*}
\begin{align*}
\int_{-\infty}^{+\infty}|\tau|^{2\gamma}|\widehat \theta^{(m)}(\tau)|^2_g
\leq & c'_6(\gamma) \int_{-\infty}^{+\infty}\frac{1}{1+|\tau|^{\alpha-2\gamma}}\|\widehat \theta^{(m)}(\tau)\|^2_{W_g}\,d\tau\\
&\qquad\qquad+c'_7(\gamma) \int_{-\infty}^{+\infty}\frac{|\tau|^{\alpha-1}}{1+|\tau|^{\alpha-2\gamma}}\|\widehat \theta^{(m)}(\tau)\|^2_{W_g}\,d\tau
\end{align*}
By Parseval inequality, the first integral is bounded as $m\to\infty$. Applying the Schwartz inequality, the second integrals yield to
\begin{align}
\int_{-\infty}^{+\infty}\frac{|\tau|^{\alpha-1}}{1+|\tau|^{\alpha-2\gamma}}\|\hat u^{(m)}(\tau)\|^2_g\,d\tau\leq&\left(\int_{-\infty}^{+\infty}\frac{d\tau}{(1+|\tau|^{\alpha-2\gamma})^2}\right)^{1/2}\\
&\qquad\qquad\times\left(\int_{-\infty}^{+\infty}|\tau|^{2\alpha-2}\|\hat u^{(m)}(\tau)\|^2_g\,d\tau\right)^{1/2}
\end{align}
\begin{align}
\int_{-\infty}^{+\infty}\frac{|\tau|^{\alpha-1}}{1+|\tau|^{\alpha-2\gamma}}\|\widehat \theta^{(m)}(\tau)\|^2_g\,d\tau\leq&\left(\int_{-\infty}^{+\infty}\frac{d\tau}{(1+|\tau|^{\alpha-2\gamma})^2}\right)^{1/2}\\
&\qquad\qquad\times\left(\int_{-\infty}^{+\infty}|\tau|^{2\alpha-2}\|\widehat \theta^{(m)}(\tau)\|^2_g\,d\tau\right)^{1/2}
\end{align}
The first integrals are finite due to $\gamma<\alpha/4$. On the other hand, it follows from the Parseval equality that
\begin{align*}
\int_{-\infty}^{+\infty}|\tau|^{2\alpha-2}\|\hat u^{(m)}(\tau)\|^2_g\,d\tau&=\int_{-\infty}^{+\infty}\|\,_{-\infty}\mathrm I_t^{1-\alpha}\tilde u^{(m)}(t)\|_g^2\,dt\\
&=\int_0^T\|\,_0\mathrm I_t^{1-\alpha} u^{(m)}(t)\|^2_g\,dt\\
&\leq\left(\frac{T^{1-\alpha}}{\Gamma(2-\alpha)}\right)^2\int_0^T\|u^{(m)}(t)\|_{V_g}^2\,dt
\end{align*}
\begin{align*}
\int_{-\infty}^{+\infty}|\tau|^{2\alpha-2}\|\widehat \theta^{(m)}(\tau)\|^2_g\,d\tau\leq\left(\frac{T^{1-\alpha}}{\Gamma(2-\alpha)}\right)^2\int_0^T\|\theta^{(m)}(t)\|_{W_g}^2\,dt
\end{align*}

Which implies that \eqref{eqgamma1} and \eqref{eqgamma2} holds. We know that subsequence of $\{u^{(m)}\}_m$ and $\{\theta^{(m)}\}_m$ (which we will denote with the same symbol) converge to some $u$ weakly in $L^2(0,T;V_g)$ and weak-star in $L^\infty(0,T;H_g)$ with $u\in L^2(0,T;V_g)\cap L^\infty(0,T;H_g)$. Similarly there exists a subsequence of $\{\theta^{(m)}\}_m$ (which we will denote with the same symbol) converge to some $\theta$ weakly in $L^2(0,T;W_g)$ and weak-star in $L^\infty(0,T;L^2(\Omega,g))$ with $\theta\in L^2(0,T;W_g)\cap L^\infty(0,T;L^2(\Omega,g))$. As $W^\gamma(0,T,V_g;H_g)$ is compactly embedded in $L^2(0,T; H_g)$ and $W^{\gamma}(\mathbb{R}, W_g, L^2(\Omega,g))$ in $L^2(0, T ; L^2(\Omega, g))$ then $\{u^{(m)}\}_m$ strongly converges in $L^2(0,T;H_g)$ and  $\{\theta^{(m)}\}_m$ in $L^2(0, T ; L^2(\Omega, g))$ respectively.

In order to pass to the limit, we consider the scalar functions $\Psi_1(t)$ and $\Psi_2(t)$ continuously differentiable on $[0, T ]$ and such that
$\Psi_1(T)=0$ ve $\Psi_2(T)=0$. We multiply (\ref{3.4})
and (\ref{3.5}) by $\Psi_1(t)$ and $\Psi_2(t)$ respectively and then integrate by parts. This leads to the equations
\begin{align*}
\int_0^T&(u^{(m)}(t),\mathrm D_{t,T}^\alpha\Psi_1(t)u_k)_g\,dt+\int_0^T b_g(u^{(m)}(t),u^{(m)}(t),\Psi_1 u_k)\,dt\\
+&\nu\int_0^T((u^{(m)}(t),\Psi_1 u_k))_g+\nu\int_0^Tb_g(\frac{\nabla g}{g},u^{(m)}(t),\Psi_1 u_k)\,dt=(u_{0m},\mathrm I_{0,T}^{1-\alpha}\Psi_2(t) u_k)_g\\
+&\int^T_0 (\xi\theta^{(m)}(t),\Psi_1u_k)_gdt+\int_0^T(f_1(t),u_k)_g\,dt
\end{align*}

\begin{align*}
\int_0^T&(\theta^{(m)}(t),\mathrm D_{t,T}^\alpha\Psi_2(t)\theta_k)_g\,dt+\int_0^T \widetilde b_g(u^{(m)}(t),\theta^{(m)}(t),\Psi_2 \theta_k)\,dt\\
+&\kappa\int_0^T((\theta^{(m)}(t),\Psi_2 \theta_k))_g\,dt+\kappa\int_0^T\widetilde b_g(\frac{\nabla g}{g},\theta_k,\Psi_2 \theta^{(m)}(t))\,dt=(\theta_{0m},\mathrm I_{0,T}^{1-\alpha}\Psi_2(t) \theta_k)_g\\
+&\int_0^T(f_2(t),\Psi_2\theta_k)_g\,dt
\end{align*}
Following the technique given in \cite{G20,T84}, as $m \rightarrow \infty$ we obtain
\begin{align}\label{eqproof1}
\int_0^T&(u(t),\mathrm D_{t,T}^\alpha\Psi_1(t)u_k)_g\,dt+\int_0^T b_g(u(t),u(t),\Psi_1 u_k)\,dt+\nu\int_0^T((u(t),\Psi_1 u_k))_g\\
&+\nu\int_0^Tb_g(\frac{\nabla g}{g},u(t),\Psi_1 u_k)\,dt=(u_{0},\,\mathrm I_{0,T}^{1-\alpha}\Psi_1 u_k)_g+\int^T_0 (\xi\theta(t),\Psi_1v)_gdt+\int_0^T(f_1(t),u_k)_g\,dt
\end{align}

\begin{align}\label{eqproof2}
\int_0^T&(\theta(t),\mathrm D_{t,T}^\alpha\Psi_2(t)\theta_k)_g\,dt+\int_0^T \widetilde b_g(u(t),\theta(t),\Psi_2 \theta_k)\,dt\\\nonumber
+&\kappa\int_0^T((\theta(t),\Psi_2 \theta_k))_g\,dt+\kappa\int_0^T\widetilde b_g(\frac{\nabla g}{g},\theta_k,\Psi_2 \theta(t))\,dt=(\theta_{0},\mathrm I_{0,T}^{1-\alpha}\Psi_2(t) \theta_k)_g\\\nonumber
+&\int_0^T(f_2(t),\Psi_2\theta_k)_g\,dt
\end{align}
This equations hold for $v$ and $\tau$ which are finite linear combination of $u_k$ and $\theta_k$, $k=1,\dots,m$ and by continuity it holds for any $v$ in $V_g$ and $\tau\in H_g$. It then follows that $\{u,\theta\}$ satisfies the two first equations \eqref{benard}. To end the proof it still to check that $\{u,\theta\}$ satisfies the initial conditions $u(0)=u_0$ and $\theta(0)=\theta_0$. To do so it suffices to multiply the two first equations in \eqref{benard} by $\Psi_1$ and \eqref{} by $\Psi_2$ respectively and then integrate. By making use of the integration by parts and comparing with \eqref{eqproof1} and \eqref{eqproof2}, one can find that
\[
(u_{0}-u(0), v)_g\,\mathrm I_{0,T}^{1-\alpha}\Psi_2(t)=0,\quad\text{and}\quad (\theta_{0}-\theta(0), \tau)_g\,\mathrm I_{0,T}^{1-\alpha}\Psi_2(t)=0
\]which lead to the desired result by taking a particular choice of $\Psi_1$ and $\psi_2$.

For the uniqueness of the weak solutions let $(u_1,\theta_1)$ and $(u_2,\theta_2)$ be two weak solutions
with the same initial condition. Let $w = u_1- u_2$ and $\widetilde{w}~ = \theta_1 -\theta_2$. Then we have

\[
\begin{array}{r}
D^{\alpha}_{t}(w,v)_g+ b_g(u_1,u_1,v)- b_g(u_2,u_2,v)+\nu (\nabla w,\nabla v  )_g+\nu (C_gw,v)_g= (\xi  \widetilde{w},v)_g
\\  \\ D^{\alpha}_{t} ( \widetilde{w},\tau)_g+ \widetilde{b}_g(u_1,\theta_1,\tau)-\widetilde{b}_g(u_2,\theta_2,\tau)    +\kappa(  \nabla\widetilde{w},\nabla\tau )_g+\kappa \widetilde{b}_g
( \displaystyle \frac{\nabla g}{g},\tau,\widetilde{w})=0
\end{array}
\]
Taking $v = w(t)$ and $\tau = \widetilde{ w}(t)$ one obtains,
\[  \displaystyle D^{\alpha}_{t}(w,w)_g+b_g(w,u_2,w)+\nu | A^{1/2}_g w |^2_g+\nu (C_g w,w)_g=(\xi\widetilde{w},w)g  \]
\[  \displaystyle D^{\alpha}_{t} ( \widetilde{w},\widetilde w)_g+\widetilde{b}_g(u_1,\theta_1,\widetilde{w})-\widetilde{b}_g(u_2,\theta_2,\widetilde{w})
+
\kappa| \widetilde{A}^{1/2}_g \widetilde{w} |^2_g+ \kappa \widetilde{b}_g ( \displaystyle \frac{\nabla g}{g},\widetilde{w},\widetilde{w} )=0.\]
By applying the bounds on the terms $b_g$,  $\widetilde{b}_g$ it then follows by Cauchy-Schwarz inequality and Gronwall like inequality that $w(t) = 0$ and $\widetilde{w}(t) = 0$ for all $t \geq 0$ since we have $w(0) = 0$
and $\widetilde{w}~(0) = 0$. Thus the Theorem is proved.

\end{proof}