International Journal Greenhouse Gas Control CO 2 storage well rate optimisation in the Forties sandstone of the Forties and Nelson reservoirs using evolutionary algorithms and upscaled geological models

Optimisation is particularly important in the case of CO 2 storage in saline aquifers, where there are var- ious operational objectives to be achieved. The storage operation design process must also take various uncertainties into account, which result in adding computational overheads to the optimisation calcu- lations. To circumvent this problem upscaled models with which computations are orders of magnitude less time-consuming can be used. Nevertheless, a grid resolution, which does not compromise the accu- racy, reliability and robustness of the optimisation in an upscaled model must be carefully determined. In this study, a 3D geological model based on the Forties and Nelson hydrocarbon ﬁelds and the adjacent saline aquifer, is built to examine the use of coarse grid resolutions to design an optimal CO 2 storage solution. The optimisation problem is to ﬁnd optimal allocation of total CO 2 injection rate between existing wells. A simulation template of an area encompassing proximal-type reservoirs of the Forties-Montrose High is considered. The detailed geological model construction leads to computationally intensive sim- ulations for CO 2 storage design, so that upscaling is rendered unavoidable. Therefore, an optimal grid resolution that successfully trades accuracy against computational run-time is sought after through a thorough analysis of the optimisation results for different resolution grids. The analysis is based on a back-substitution of the optimisation solutions obtained from coarse-scale models into the ﬁne-scale model, and comparison between these back-substitution models and direct use of ﬁne-scale model to conduct optimisation. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Geological storage of CO 2 is regarded as a viable option for tackling climate change. The main types of geological systems considered for CO 2 storage are depleted oil and gas fields, saline aquifers and unmineable coal seams. In order to achieve long term trapping of the CO 2 in the storage system, it is necessary that an impermeable caprock exists above the storage formation and the geological structure allows for the CO 2 to be confined within the intended formation. In the actual design of CO 2 storage site operation it is important that potential risks are minimised or entirely avoided.
One such potential risk is CO 2 migration outside the targeted geological storage complex. This might occur if the CO 2 plume is allowed to reach the structural spill point of the targeted zone. Moreover, the CO 2 injection should not affect regions outside the licensed area boundaries and there should be no migration of CO 2 into the neighbouring fields. The crucial objectives of CO 2 storage are to utilise the storage capacity optimally, contain the CO 2 as a trapped phase and avoid the risks of leakage or spillage from the storage complex. These competing requirements necessitate that optimisation techniques are used in evaluating the performance of the proposed CO 2 injection system designs.
The objective of this paper is to illustrate the use of optimisation techniques in conjunction with upscaled models of heterogeneous storage complexes in the North Sea, with the aim to minimise the risk of lateral migration of injected CO 2 reaching the system spill point, while maximising the amount of CO 2 stored. One important aspect considered in this process is that the upscaled static models should account for reservoir heterogeneity while ensuring that http://dx.doi.org/10.1016/j.ijggc.2016.04.011 1750-5836/© 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). the detailed system is represented adequately. The storage complex used to illustrate the approach comprises reservoirs of the Forties and Nelson oilfields and their surrounding aquifer within the UK Central North Sea. These two fields have already been cited in the literature for high potential for CO 2 storage (Cawley et al., 2005;Ketzer et al., 2005;SCCS, 2009) or CO 2 -EOR (Kemp and Kasim, 2012). The two fields are assumed to be fully depleted and saturated with brine from the surrounding aquifers. The optimisation study presented here assumes that a fixed daily supply rate of CO 2 is available for injection. The optimisation variables are the injection rate at each well, namely the allocation of rate to each well from the total available CO 2 . The optimisation objective is to minimise the fraction of CO 2 that is in free gaseous state outside the licensed regions and to maximise the fraction of CO 2 that is stored as dissolved CO 2 in the aqueous phase or residually trapped CO 2 in the gaseous phase.
An evolutionary optimisation algorithm is chosen to conduct the work as the objective functions are multimodal and the calculation of gradient or the use of a gradient-based optimisation algorithm would not be appropriate. Since the geological modelling of the site is very detailed and computationally expensive for CO 2 storage site assessment, the study focuses on the interplay of single and multiobjective Genetic Algorithm (GA) optimisation and upscaling for CO 2 storage. The upscaling method used is a numerical single-phase permeability upscaling. The optimisation result of upscaled model is exercised on the fine-scale grid, in a back-substitution procedure, to determine the reliability of the upscaled models for use in the optimisation. The outputs of optimisation from three upscaled levels of the realistic channelised aquifer model are compared with the outputs of the original fine-scale model and the effect of averaging of sub-grid flow on the optimisation output of the coarse-scale grid is assessed.

Background
Numerical simulation of CO 2 storage processes has been used widely to estimate the storage capacities of various sites (Doughty and Pruess, 2004;Flett et al., 2007;Kopp et al., 2009 andZhang, 2008). The estimation of CO 2 storage efficiency, i.e. the fraction of the available pore space that is utilised for CO 2 storage (van der Meer, 1995), is a crucial step in the evaluation of a candidate CO 2 storage site. Moreover, to achieve the highest storage efficiency and containment of CO 2 within the geological time-spans, it is important to fully understand and assess the mechanisms of CO 2 storage, which have been studied in great detail by Gunter et al. (1997), Holtz (2002, Xu et al. (2003), Holloway (2005) andIPCC (2005).
The dynamic methods, such as those developed by Doughty and Pruess (2004), Flett et al. (2007), Kopp et al. (2009) and Pruess and Zhang (2008), involve numerical simulation of physical phenomena that include a number of different mechanisms of trapping. These methods are more accurate than static methods and are widely applicable to real cases as they account for various flow dynamics, including multi-phase flow, density-dependent flow, heat transfer, ground-water hydrodynamics, CO 2 dissolution kinetics and geochemical reactions, which are fully described by the finite-element or finite-difference discretisation of the partial differential equations. More importantly, numerical models take into account the heterogeneity of the subsurface structures which is an important parameter influencing both migration patterns and trapping mechanisms for CO 2 injected in open saline aquifers (Gasda et al., 2013).
The detailed discretisation in large simulation models, however, leads to large number of unknowns and the computations become intensive (Aarnes et al., 2007;Matthäi and Nick, 2009). Furthermore, to account for the uncertainty in model parameters − such as permeability, porosity, temperature, brine salinity, fracture properties and relative permeability − multiple realisations should be considered. This adds significantly to the computational burden of CO 2 injection design so as to utilise the storage capacity. One well suited solution is upscaling. Upscaling has been used in the context of CO 2 storage in a number of studies such as Behzadi and Alvarado (2012); Gasda et al. (2013); Green and Ennis-King (2010); Juanes and MacMinn (2008); Mouche et al. (2010) and Saadatpoor et al. (2011).
However, upward migration of CO 2 has been particularly challenging in upscaling. Saadatpoor et al. (2011) studied the effect of permeability upscaling on buoyancy-driven vertical flow of CO 2 and the mechanism of CO 2 immobilisation by local trapping. Different degrees of coarsening were considered and, based on average gas saturation and mass of CO 2 in the storage aquifer, the upscaling methods were found to underestimate the extent of CO 2 trapping because upscaling averaged the map of capillary pressure. Green and Ennis-King (2010) presented simple analytical expressions for the mean and variance of the upscaled vertical permeability distribution in a reservoir with randomly distributed impermeable barriers. Comparing the heterogeneous reservoir and a homogeneous reservoir with effective vertical permeability, CO 2 reached the top of the formation sooner for the heterogeneous medium. Additionally, convection began much sooner in the heterogeneous cases than in homogeneous case with the same effective permeabilities. A number of studies, such as that of Gasda et al. (2011), Juanes and MacMinn (2008) and Mouche et al. (2010) also focused on vertical-direction upscaling for CO 2 storage by vertical equilibrium assumption. Given the strong buoyancy forces, it is reasonable to assume that complete gravity segregation occurs quickly during and after the injection period. In addition, the large horizontal and thin vertical scales result in instantly developing vertical movement of the fluids (Gray et al., 2012), as if the vertical permeability is actually infinite. This assumption facilitates vertical integration of the governing three-dimensional flow equations to obtain a set of two-dimensional equations.
When the viscous force dominates the buoyancy, the upscaling methods that are used for oil reservoirs are considered to be applicable for upscaling CO 2 storage models in the lateral direction. For instance, Li et al. (2011) used a single-phase steady state permeability upscaling method to create a coarse-scale model. The objective was to derive a resolution at which the coarse-scale model is sufficiently accurate for predicting gas migration, CO 2 storage, and pressure build-up in the absence of capillarity and viscousdominated flow.
The difference between the comparative and evaluative analyses of upscaling for aquifer models and oil and gas reservoirs is that, for the latter, hydrocarbon production curves throughout the simulation exist and can be employed as a criterion of accuracy of the upscaled model. Simply put, the closer the curve obtained by the upscaled model is to that of the fine scale model, the more reliable the upscaled model is. However, this is not the case for the aquifers with no hydrocarbon production data. Therefore, other evaluative criteria should be employed. One criterion suggested and employed in this study is the performance of the upscaled models in optimisation. The convergence curves of the optimised solutions, the Pareto fronts of the multiple objective optimisation and the solutions are some of the outputs of the optimisation that are used for evaluation of the upscaled models in this work.
A number of researchers (Kovscek and Cakici, 2005;Bergmo et al., 2011;Leach and Mason, 2011;Cameron and Durlofsky, 2012;Shamshiri and Jafarpour, 2012;Babaei et al., 2015) have focused on single and multi-objective optimisation for CO 2 storage with different optimisation variables and performance measures. As the computationally expensive design procedure involved in optimising large subsurface models is a major concern, the use of upscaled models is more justified. In this line of research, Cameron and Durlofsky (2012) used a modest level of coarsening (by a factor of two in each direction) in their non-gradient-based optimisation and assessed the impact of grid resolution on the optimised solutions. They observed a reasonable correspondence in mobile CO 2 fraction between the models. Moreover, a study of gridding effects on CO 2 storage simulation has been performed by Yamamoto and Doughty (2011). The conclusion was that there is an underestimation of gravity override and maximum plume extent by a coarse grid. They also recommended near well upscaling for capturing salt precipitation near the wells.
The upscaling of the petrophysical model in this work is limited only to the lateral direction. One reason for not upscaling vertically is that the vertical resolution is more sensitive to upscaling in estimating the upward migration of CO 2 . The other reason is that we only use one thin (1 meter thickness) layer of the geological model, representing a geological zone in a large multi-zonal aquifer with pressure and flow discontinuities in the vertical direction between the zones. The objective functions are dependent on the extent of CO 2 plume and the CO 2 storage mechanisms. The static model of the Central North Sea Paleocene/Eocene Forties Sandstone Member used in this study provides a sufficiently complex environment for evaluation of the upscaled models.

Description of the simulator and upscaling methodology
The compositional simulator ECLIPSE E300 Schlumberger, 2010 is used to conduct the reservoir simulations. The CO2STORE option is activated to simulate injection and geological storage of CO 2 into the saline aquifer systems that were modelled. The estimates from the simulation are the molar mass of CO 2 that is dissolved, as well as the residually trapped and free CO 2 found in the gaseous state. For the simulations using the CO2STORE option, three phases are considered, a CO 2 rich phase, a H 2 O rich phase and a solid phase. The mutual solubilities of CO 2 and H 2 O are calculated to match experimental data for CO 2 -H 2 O systems under typical CO 2 storage conditions ranging between 12 and 100 • C and up to 600 bar. These are calculated following the procedure given by Spycher and Pruess (2005), based on fugacity equilibration between water and a CO 2 -phase. Water fugacity is obtained by Henry's law, while CO 2 fugacity is calculated using a modified Redlich-Kwong equation of state (Redlich and Kwong, 1949). The components that can be considered in the CO2STORE option are CO 2 and H 2 O present in water and gas phases and the salts NaCl, CaCl 2 and CaCO 3 present in either aqueous or solid phase Schlumberger, 2010. The upscaling algorithm used for porosity and the net-to-gross ratio is the arithmetic volumetric averaging. For permeability, the single-phase non-tensorial pressure solver method (Begg et al., 1989;Christie, 1996) is used, which is a simple method. More advanced upscaling methods for highly heterogeneous permeability fields such as methodologies presented by Chen and Durlofsky (2006) and Evazi and Jessen (2014) could be implemented. The main expression for the upscaled permeability or the block equivalent permeability of coarse gridblock E in the pressure solver method, denoted henceforth by K * E , is given by: where p is the pressure and V represents the volume of gridblock E. This equation can either be solved numerically or, in special cases, it has an exact solution. The pressure solver method aims at obtaining and inverting the solution of Eq. (1) (calculated locally) over the domain of coarse gridblock E, to derive K * E . The solution is dependent on the choice of boundary condition. One common choice utilised here is to assume a generic axis-oriented boundary condition (two sides with a prescribed pressure gradient and two sides with no flow). Then, if x is the direction of the pressure gradient, the total flow rate q c x is computed by summation of outlet fluxes, thus, the diagonal element of where p is the assumed pressure gradient, x is the length of the gridblock in x direction, A is the area from which the outlet fluxes exit. By alternating the boundary condition over the sides, K * E is estimated for the y and z directions.

Evolutionary optimisation
In general, a wide variety of optimisation techniques using gradient descent methods may be used if the objective function is smooth and analytically tractable. However, this is not commonly the case for subsurface models. Additionally, if the functional form is known to be convex, then semi definite programming or other linear-matrix-inequality-based techniques will be more expedient, as commercial solvers that solve such problems in polynomial time exist.
From the optimisation point of view, simulators such as ECLIPSE needs to be treated as black box models. Therefore, the traditional methods of optimisation are not suitable in general. Intelligent bioinspired optimisation techniques like Genetic Algorithm (GA) or Particle Swarm Optimisation are well suited to this kind of problems and, unlike the traditional gradient-based optimisers, can handle discontinuous, noisy and stochastic objective functions.
The GA encodes the solution variables as genes and initialises a population randomly within the lower and upper bounds of the search space. The fitness of each of the genes is evaluated by using the user specified objective function. The algorithm then uses crossover and mutation operators to evolve the next generation of the population. A few elite genes with the highest fitness function values are directly passed on to the next generation to preserve the best obtained solutions. The algorithm iteratively performs these operations to produce a newer population of genes until a specified number of generations are completed. The best gene in the last generation represents the optimised solution obtained by the algorithm.
The multi-objective non-dominated sorting GA algorithm (Deb et al., 2002) is used in this study to produce the Pareto front solution. The algorithm uses the same basic principles of mutation and crossover as the single objective GA. However it utilises a sorting scheme to rank the Pareto fronts in the order of their nondomination. A crowding distance assignment algorithm is used which gives more importance to solutions in sparsely populated regions of the Pareto front and, hence, helps in obtaining a better spread of the Pareto solutions. These mechanisms assist in applying selection pressure on the solutions and the approximate Pareto front evolves towards the true Pareto front over the generations.

Objective functions for optimisation
The optimisation problem is defined using two objective functions: where are the dissolved, residually trapped and mobile CO 2 present in gridblock i at time t respectively, therefore J 1 is the fraction of injected CO 2 as free gas present in the M b gridblocks of region M b ⊂ N b at the end of simulation. Region M b is the region of the aquifer lying outside the licensed regions of the Forties and Nelson reservoirs and N b is the entire gridblocks of the system. The second optimisation function, J 2 , is the fraction of injected CO 2 stored by dissolution or residual entrapment in the entire N b gridblocks of the domain, calculated as the summation over time of simulation (t max ). In the multi-objective optimisation setting, the goal is to minimise free gas outside the boundaries of the reservoir (J 1 ) and maximise trapped CO 2 (J 2 ). For the single objective case only J 1 is minimised. The reservoir simulator outputs (v i,t ) depend on the total flow rate, Q , and the allocation vector of injection rates, x, for the wells that are present in the system.
The dissolution and residual trapping are enhanced by accessibility of fresh brine to gaseous CO 2 . A better sweep efficiency helps this process and can be achieved through better engineering design of CO 2 injection. However, with larger extent of CO 2 plume and maximised sweep efficiency, leakage or spillage may induce a risk of CO 2 migrating outside the licensed storage areas. In other words, the objective functions defined are in fact contradicting each other.

The study area
The site modelled is located on the Forties-Montrose High in the UK Central North Sea and includes the Forties and Nelson hydrocarbon fields ( Fig. 1; Wills , 1991). The hydrocarbon reservoir consists of submarine fan deposits of the Paleocene/Eocene Forties Sandstone Member overlain by Lower Eocene shales (Hughes et al., 1990;Whyatt et al., 1992). The reservoir is located in the proximal inner (interbedded sand/shale) to middle (mainly massive sand) part of the Forties Fan system (Hughes et al., 1990) and is mostly channelised and characterised by high net to gross sandstone ratios, good porosities and high permeabilities (Hempton et al., 2005).
The reservoir comprising the geological model was subdivided into 8 Zones, Zone M, the overlying top seal, an 'Upper Sand' (Zones L, K, J, H, F, and E) separated by a continuous and uninterrupted mudstone layer from the 'Lower Sand' (Zone D). The mudstone layer is a barrier to fluid flow within the reservoir and results in a vertical pressure discontinuity between the Upper and Lower sands (Wills and Peattie, 1990). In addition, an extremely effective and laterally extensive permeability barrier, referred to as the 'Charlie Shale' (Denny and Heusser-Maskell, 1984) exists within the upper part of Zone H that effectively separates the Upper Sand Zones L, K, J above (only existing in the western parts of the Forties Field) and Zones H(lower part), F, and E beneath. Pressure communication within the Zone J overlying the Charlie Shale is good, with very little indication of vertical permeability barriers. In the reservoir zones beneath the Charlie Shale, minor permeability barriers have been indicated. The flanks of the structure have fairly gentle dip over the reservoir section (Hughes et al., 1990). It is assumed that the reservoirs are fully depleted and saturated with aquifer brine before CO 2 injection starts.
Reservoir pressure in both fields was initially maintained by basal aquifer influx (Simpson and Paige, 1991). Since 1976, aquifer support has been supplemented by peripheral seawater injection (Mitchell, 1978) so that the pressure was sustained throughout the oil production stage. With such evidence, the assumption that the whole system is an open boundary aquifer is appropriate.

Geological model and attribution
The geology of the Forties and Nelson hydrocarbon fields in the UK sector of the North Sea has been widely studied (e.g., Kulpecz and van Geuns, 1990;Kunka et al., 2003). The geological model developed in this study broadly captures and represents the het- four-way dip-closed anticlinal structure comprise the main hydrocarbon producing fairways of the two hydrocarbon fields (Fig. 2(a)). The amalgamated channels and their associated channel margins and interchannel areas migrated laterally through time resulting in the variation in relative dominance and position of the different parts of the submarine fan system and a high degree of lateral and vertical variation in the reservoir. This is represented in the lithologies found in the system and their associated petrophysical parameters.
The geological modelling was carried out in three parts: (a) modelling the amalgamated channels in the fairways; (b) modelling the interchannel areas; and (c) attributing the model with petrophysical properties. A typical cross-section of amalgamated channels and channel margins consists of the channel sands, low permeability basal lags, high permeability basal lags and the intrachannel doggers.
Considering each zone of the model separately, the first step in building the static model was to construct a reasonable fairway structure in the zone using the interpreted channel boundaries. Fig. 2(b) shows the flow line for channels in one zone, constructed to represent the maximum sediment transport line. This was followed by modelling of the amalgamated channels inside the fairway volume using Object Modelling in the PETREL software and the attributes specified in Table 1 for each of the zones. This allowed to stochastically model the channel complex based on parameters such as facies associations, the number or percentage of channels and shale doggers, channel layout and section parameters. Fig. 2(b) illustrates the facies associations in Zone E. In order to model the interchannel slump bodies, the fairway polygons were inflated to include buffer areas around the poly-gons, assuming distances varying between 2 and 5 km. Ellipsoidal shapes were assumed typical for slump bodies owing to their generally elongated shapes. Table 2 summarises the ellipsoidal body parameters used.
The attribution of the geological model with petrophysical properties, namely porosity, permeability and net-to-gross (NTG) ratio, was carried out using Gaussian random functions (Lifshits, 1995). The ranges of values used, including their mean values, are summarised in Table 3 for the different geological facies. Some of these values are based on generalisations from sparse literature references such as Kunka et al. (2003) and Wills (1991) and, in the absence of data, some assumptions were made by the authors based on geological expert judgement.

Upscaling
The geologically detailed fine-scale model is laid out on a 844 × 640 lateral resolution grid with an average gridblock length of 50 m in x and y directions. An average thickness of one metre is used for layering the model. The geological model is referred to as the Fine-scale Model (FM) from this point on. This resolution is appropriate as it allows us to describe the channel and slump deposits' extent and the lateral and vertical heterogeneity within channels. As vertical heterogeneity within the individual channels is much more significant, considering the mode of presence for levee and shale dogger facies, it was decided to use the maximum of  Table 3 Petrophysical properties for different facies types from Kunka et al. (2003) and Wills (1991)  width variation for basal lags and doggers and much finer thickness for the individual facies objects.
Here only the top layer of Zone E that comprises 392,155 active gridblocks is chosen for optimisation study. This layer is upscaled in x and y directions into three coarser resolution girds to have 4 resolution grids in total: As mentioned the directional absolute permeabilities (Permeability I, J, K in x, y and z directions respectively) are upscaled by flow-based upscaling whereas porosity and net-to-gross ratio (NTG) are upscaled by volume averaging. To illustrate, Permeability I and Permeability J for FM and CM3 are shown in Fig. 3, and Porosity and NTG for FM and CM3 are shown in Fig. 4. The upscaling has made no visually identifiable difference to the petrophysical properties of interest.

Injection simulations
The eight wells used for CO 2 injection simulation are shown in Fig. 2(b). Two of the four wells in Forties (vertical wells) were used for peripheral water injection during oil production and two others (deviated wells) were used to produce oil at platform FC (Forties Charlie) located in the western flank of the oilfield. Four wells in Nelson (all deviated) were all oil production wells drilled from a single platform located in the middle of the oilfield. These wells were selected because they were drilled in the channel area of Zone E.
The brine properties used in simulations assume an initial molar composition of 0.05 NaCl, 0.04 CaCl 2 and 0.91 water. The normal diffusion coefficients were set to 0.001 cm 2 /sec for water and CO 2 components in the gaseous phase, and 0.0001 cm 2 /s for all components present in the aqueous phase. CO 2 density at gaseous state is obtained by Redlich-Kwong equation of state which is tuned following Spycher and Pruess (2005) in ECLIPSE E300 to accurately give the density of the compressed gas phase. The CO 2 gas viscosity is calculated in ECLIPSE E300 from Vesovic et al. (1990) and Fenghour et al. (1998).
Formation brine properties, either obtained from literature or in the absence of data set to ECLIPSE E300 default values, are summarised in Table 4. It should be noted that the properties that exhibit different values for each oilfield are averaged and applied for the whole aquifer study area.
The relative permeability curves, shown in Fig. 5, were derived for gas-water system from a three phase oil-brine-CO 2 system in an investigation of CO 2-EOR potentials for a segment of the Forties field by Cawley et al. (2005). Based on these curves, connate water saturation (S wc ) is 0.2 but the relative permeability of water does not increase significantly until water saturation reaches 0.4. This implies that CO 2 saturation can achieve large values of maximum gas saturation (Sg max = 1 − S wc = 0.8) only near wells. In the absence of field data, the capillary pressure is neglected and is set to zero and no hysteresis is considered for the relative permeability curves.
The CO 2 storage license regions were assumed to be the reservoirs. It should be noted that the reservoir boundaries used are approximate and are not considered accurate delineations of either depth or fluid contacts. Although the particular choice of boundaries affects the optimal solution obtained as a result of the analysis,

Table 4
Formation brine properties.

Property
Forties Nelson Salinity (ppm of NaCl) 55,500 (Wills, 1991) 84,000 (Kunka et al., 2003) Resistivity (ohm m) 0.034 (Wills, 1991) 0.034 Formation volume factor (rm 3 /sm 3 ) 1.24-1.32 (Wills, 1991) 1.357 (Kunka et al., 2003 C) 96 (Wills, 1991) 107 (Kunka et al., 2003) Initial fluid pressure (bar) 222 (Wills, 1991) 229 (Kunka et al., 2003) it does not affect the applicability of the combined modelling and optimisation methods proposed in this paper. A total rate of 0.7 million sm 3 per day, approximately equivalent to 0.5 million tonnes of CO 2 per year is injected into the aquifer for 30 years and an additional 50 years are allowed for stabilisation and investigation of possible migration of CO 2 . During the injection period, a total of 15 million tonnes of CO 2 are injected.
Although, the time period investigated for the aquifer injection and stabilisation in the geological sense is limited, nonetheless, this CO 2 storage model is adequate to study the trade-offs in the design and the effect of upscaling in the optimisation process, which is the objective of the methodology illustrated here. Larger scale models and simulation over longer injection periods would only require additional simulation time.
The simulations were run on a Windows based operating system with an Intel Core i7-4820K 3.70 GHz CPU and 32.0 GB RAM. Each of FM, CM1, CM2 and CM3 takes 4830, 555, 179 and 96 s to run, respectively. These numbers must be multiplied by 800 (20 populations × 40 generations) to indicate the time required for optimisation of each model and each grid. As shown in Fig. 6, the difference in run times between FM and CM1 is higher than run times between CM1 and CM2 and between CM2 and CM3, suggesting that the run time does not scale linearly with the size of the model. As the optimisation for FM requires 800 × 4830 s or in other words 44.72 days, we used parallel computing allowing the runtime reduction to order of 6 days instead. However the parallelisation has not taken into account in Fig. 6.

Single objective optimisation
Starting with the single objective optimisation, the convergence profiles of the final optimum solution from one generation to    another for different grids are shown in Fig. 7(a). The convergence is assured only when the best fits do not vary significantly over the iterations. It is obvious that CM3 with the coarsest resolution grid has led to loss of variation in the solution for objective function and consequently the profile is flat with no comparability of the final solution after 40 generations with that of FM, CM1 and CM2. In contrast, CM1 and CM2 have preserved the variation and "optimisability" of the objective function compared to FM and as such the final values for J 1 after 40 generations is close to that obtained by FM.
It is interesting to see the change in the order compared to FM increases from CM1 to CM3, with CM3 only having W6 as the highest share-taker similar to FM. In order to assess the degree of reliability of the coarse-scale solutions of optimisation obtained from CM1, CM2 and CM3, these solutions are used in fine-scale simulation. In other words we substitute the solutions onto FM resolution grid. These substituting simulation are denoted as CM1in-FM, CM2-in-FM and CM3-in-FM hereafter. Fig. 8 shows how optimised solutions reduce the fraction of mobile CO 2 outside the licensed regions (J 1 ) at the end of simulation. Compared to the base case scenario-where the total rate is equally divided between the eight wells, the final result of optimisation, referred to as optimal solutions, reduce fraction of mobile CO 2 from 0.2733 to 0.2166 for FM. In terms of convergence of the results from coarse resolution grids to fine resolution grid, J 1 at the end of simulation for CM3, CM2, CM1 and FM are respectively 0.2871, 0.2710, 0.2764 and 0.2733, showing a relative error of 5% for CM3 to around 1% for CM1 and CM2. Also shown in Fig. 8 are the values of J 1 throughout the simulation obtained by CM1-in-FM, CM2-in-FM and CM3-in-FM cases. The results are very reassuring, showing that the CM1-in-FM and CM2-in-FM cases both are "almost" reproducing the FM optimal solution: J 1 from CM1-in-FM at the end of simulation is 0.2224 and J 1 from CM2-in-FM at the end of simulation 0.2226. This confirms that the solutions for the injection rates shown in Fig. 7 for CM2 and CM3 both are actually optimal for FM. The solution for CM2 is of course less computationally intensive to obtain than that of CM1. Finally J 1 from CM3-in-FM at the end of simulation is as high as 0.2384, and therefore the solution can be deemed sub-optimal.
Finally, Fig. 9 illustrates the CO 2 saturation for base case and the optimal solution at the end of injection time that exhibits the larger spread of mobile CO 2 outside the license area for base case northern and southern boundaries of Forties. In the optimal case a higher injection share for W6 has led to a better filling-up of Nelson compared to base case scenario.

Multi-objective optimisation
For the multi-objective optimisation, the pareto fronts of different resolution grids, with 20 populations per generation and 40 generations, are shown in Fig. 10(a). In addition to four sets of Pareto fronts and feasible (sub-optimal) solutions for FM, CM1, CM2 and CM3, the substitution of CM1, CM2 and CM3 pareto front solutions onto FM resolution grid are conducted, and denoted by CM1-in-FM, CM2-in-FM and CM3-in-FM cases. The results are compared with FM pareto front in Fig. 10(b).
The figure shows: • The conflicting situation at which maximising the fraction of dissolved CO 2 (J 2 ) will lead to an unfavourable increase in the fraction of mobile CO 2 outside the licensed region (J 1 ). • The coarse resolution grids lead to a systematic overestimation of J 2 , and as such the dissolved amount of CO 2 will be more and more unreliable by increasing the coarsening factor. • The optimisation problem defined in this paper has not allowed a significant variation in J 2 (around 1 percent difference between lower and upper bounds of the pareto fronts in any of the resolution grids.) • The substitution cases in Fig. 10(b) show sub-optimal situations for CM3-in-FM, whereas the substitutions from CM1 and CM2 are regenerating the pareto front of FM by negligible error. Unfortu-nately there is no way to define an error metric between these points as they are not compared one-to-one and the number of pareto front solution vary from one resolution grid to another. • As of single objective optimisation, CM2 is showing a good level of reliability compared to CM1, and since the optimisation with CM2 entails less computation intensity, this resolution grid is more favourable.

Conclusions and future works
This study showed that the use of upscaled models for the optimisation of well rate allocation in designing CO 2 storage operations is feasible only if a careful preliminary assessment of the upscaled models' performances against that of geologically detailed fine scale model is carried out. A methodology was proposed for determining the accuracy of upscaled models based on back-substitution of the solutions of the optimisation from the upscaled model into the fine grid model. This procedure was used to assess the reliability of the upscaled models for the optimisation and design of a CO 2 storage operation in a heterogeneous aquifer system in the North Sea. Our results for both single and multiple objective optimisation showed that a resolution grid of 150 m by 150 m could be used safely to run optimisation reliably replacing the geologically fine resolution of 50 m by 50 m. This conclusion is only specific with the facies types of the present model and could not be readily extended to other geological models or even other geometrical systems with different vertical resolutions.
There are a few issues that remain to be addressed and can be pursued in future research. We only considered a two-dimensional model and therefore the effects of vertical resolution in upscaling were not considered. Uncertainty has not been thoroughly investigated in the present study. A methodology to handle uncertainty in the optimisation framework needs to be considered and the effect of these upscaled models with the quantification of uncertainty needs to be investigated. Other surrogate modelling techniquessimilar to the use of upscaled models -to reduce the run time complexity of these large aquifer models can be studied and compared with the upscaled models. There are several other performance objectives that may be considered for the design of the storage site. For example, from the licensing point of view, minimising the pressure build-up can be a potential objective in optimisation. Moreover, CO 2 -EOR, where oil production can be used as another objective function, is an attractive field of study to investigate the applicability of the upscaled models.
Finally, the state variable for the optimisation used in this work was the allocation of a fixed total CO 2 injection rate. In a complete operational design of a storage site many other design parameters, such as the number and locations of wells, the maximum allowable bottomhole pressure during injection, need to be defined. In all such cases, the upscaling needs to be quantitatively re-evaluated and the approach presented in this work can be used for this purpose.