Optimal design and operation of integrated wind-hydrogen-electricity networks for decarbonising the domestic transport sector in Great Britain

Abstract This paper presents the optimal design and operation of integrated wind-hydrogen-electricity networks using the general mixed integer linear programming energy network model, STeMES (Samsatli and Samsatli, 2015). The network comprises: wind turbines; electrolysers, fuel cells, compressors and expanders; pressurised vessels and underground storage for hydrogen storage; hydrogen pipelines and electricity overhead/underground transmission lines; and fuelling stations and distribution pipelines. The spatial distribution and temporal variability of energy demands and wind availability were considered in detail in the model. The suitable sites for wind turbines were identified using GIS, by applying a total of 10 technical and environmental constraints (buffer distances from urban areas, rivers, roads, airports, woodland and so on), and used to determine the maximum number of new wind turbines that can be installed in each zone. The objective is the minimisation of the total cost of the network, subject to satisfying all of the demands of the domestic transport sector in Great Britain. The model simultaneously determines the optimal number, size and location of each technology, whether to transmit the energy as electricity or hydrogen, the structure of the transmission network, the hourly operation of each technology and so on. The cost of distribution was estimated from the number of fuelling stations and length of the distribution pipelines, which were determined from the demand density at the 1 km level. Results indicate that all of Britain's domestic transport demand can be met by on-shore wind through appropriately designed and operated hydrogen-electricity networks. Within the set of technologies considered, the optimal solution is: to build a hydrogen pipeline network in the south of England and Wales; to supply the Midlands and Greater London with hydrogen from the pipeline network alone; to use Humbly Grove underground storage for seasonal storage and pressurised vessels at different locations for hourly balancing as well as seasonal storage; for Northern Wales, Northern England and Scotland to be self-sufficient, generating and storing all of the hydrogen locally. These results may change with the inclusion of more technologies, such as electricity storage and electric vehicles.


Introduction
The advantages of hydrogen as an environmentally clean fuel can be fully realised when it is produced from renewable energy sources. Hydrogen is a flexible energy storage medium that can be used for both short- 35 and long-term storage applications, in addition to being a versatile intermediate that can be converted to electricity, heat and transport fuel. Hydrogen, like electricity, can complement renewable sources particularly well -both are energy carriers that can transmit energy from primary energy sources to end-users.
It is widely accepted that hydrogen may have a role to play in decarbonising the transport sector, which still relies almost exclusively on oil. In Great Britain (GB), for example, the domestic transport sector is a 40 major oil user and is responsible for approximately 20% of total GB carbon dioxide emissions [2]. Indeed, decarbonising this sector is a main driver behind the development of fuel cell and electric vehicles. On a more positive note, GB has a very strong potential for wind power; in fact, it is considered one of the best locations in the world and the best in Europe [3]. Converting wind energy to either electricity or hydrogen that can be used in electric or fuel cell cars results in zero emissions (or low emissions if the emissions in 45 manufacturing and installing the network components are considered), which can help to achieve future emissions targets [4,5,6].
Since both demands and wind availability are distributed in space and vary with time, there is no guarantee that wind power will be available where and when it is needed. Therefore, a network of technologies that is sufficiently flexible to deal with the mismatch between the intermittent supply and demand for energy 50 is needed. It is then natural to ask what the network will look like. How many wind turbines will be needed and where will they be located? What role will electricity networks have to play: should the energy generated by wind turbines be transmitted as electricity or converted in-situ to hydrogen and transmitted through pipelines? Even at the national level, where different regions are interconnected by transmission lines, there may not be sufficient wind power to meet peak demands. Therefore, storage technologies are 55 expected to play a role; their type, size and location need to be determined. Conversion technologies, such as electrolysers and fuel cells, to interconvert electricity and hydrogen may also be necessary.
The questions posed above are extremely difficult to answer without using mathematical models, because the networks can be highly complex and integrated with many alternative options. Several models for hydrogen networks have been developed and typically fall into one of the following categories: focus on design with area occupied by technologies (which for wind turbines was considered significant). Finally, a more detailed model of the hydrogen network was included, where hydrogen is generated and stored at a high pressure and transmitted and distributed at a lower pressure. The compressors and expanders required to move the hydrogen from one pressure level to another were also included in the model of the hydrogen network, along with their interaction with the electricity network (i.e. electricity consumed by the compressors and 130 electricity generated by the expanders).
To date, to the authors' knowledge, the model presented in this paper is the first MILP model in the literature for integrated wind-hydrogen-electricity networks that can simultaneously determine the design and operation of the network while considering both spatial and temporal aspects in detail so that transport and storage of energy can be modelled more accurately. 135 The paper is structured as follows: Section 2 defines the problem and Section 3 characterises the properties of the system based on the spatial and temporal discretisation used in the model. The wind farm siting analysis, performed using GIS, is also presented in Section 3. Section 4 describes the structure and components of windhydrogen-electricity networks, the mathematical model for which is presented in Section 5. The modelling of the distribution network is discussed in Section 6. The results of the case studies are presented in Section • The characteristics of each technology, e.g. efficiency and unit costs (capital, operating and maintenance costs) Determine: • The optimal number, size and location of wind turbines, electrolysers, hydrogen storage, fuel cells, • Minimise total network costs 3. Spatio-temporal representations 160 This section characterises the properties of the system according to the spatial and temporal representations used in the model.

Spatial discretisation
STeMES models the spatial dependencies of the system by dividing the study region into a number of zones, each of which is assumed to have uniform properties (e.g. resource demand and availability) and may host 165 a certain number of technologies for generation/conversion and storage. Infrastructures for transporting resources may connect each zone to each of its neighbours.
In this study, Great Britain was divided into 16 zones based on the National Grid Seven Year Statement (NG SYS) 17 study zones [33]. Figure 1 shows the 16 transmission zones considered in this study: Z1 to Z3 correspond to the same zones in the NG SYS; NG SYS's zones 4 and 5 were combined to form Z4 because 170 zone 5 is much smaller than the other zones and keeping the number of zones to a minimum helps to reduce the computational burden of the model; and Z5 to Z16 correspond to NG SYS's zones 6 to 17.
One advantage of using a similar spatial discretisation to that of the National Grid is that a significant amount of data, in both zonal and national form, are available from sources such as [34]. Where more detailed data are available, they are aggregated for each zone.

Wind turbine siting constraints
In each zone, only a certain amount of land is available for siting of wind turbines. The maximum land area in each zone was determined by applying a number of technical and environmental constraints. The technical constraints take into account the site's wind speed, topography and accessibility, whereas the environmental constraints consider the landscape impacts and planning restrictions [35,36,37]. In this work, the following 180 10 criteria were used to determine the total land area in each zone suitable for siting wind turbines using GIS.
1. Wind speed An average annual wind speed of at least 5m/s measured 45m above ground level is needed to justify the installation of wind turbines on economic grounds [35,36,37]. 185 2. Slope Sites with slope of less than 15 percent were considered suitable for wind turbines to ensure that the parts can be safely transported for installation on steep mountainous areas [38].

Access
A maximum distance of 500m from the minor road network is imposed to allow entry of construction 190 vehicles, delivery of materials and general access for supplies and staff [35,36].

Connectivity to National Grid
The wind turbines need to connect to an energy distribution/transmission network. For simplicity, it was assumed that National Grid lines closely follow the road and that new distribution/transmission lines, if they are to be built, will be laid along existing lines. Therefore, a suitable site should not be 195 more than 1500m from the main road [35,36]. For safety, a buffer of 200m from the main road was also applied.

Planning restrictions
Locations that are nationally designated as nature and science protection areas were excluded. This was done by excluding the areas classified as Sites of Special Scientific Interest (SSSI) [35,36]. 200

Population impacts
For safety and to minimise noise intrusion, a buffer of 500m from developed land used area (DLUA) was applied [35,36].

Water pollution
To minimise the impact on wildlife in and around water courses and to decrease the risk of flooding 205 of the wind turbine during winter and spring, only sites that are more than 200m from a river were considered suitable [35,36].

Interference
To minimise impact on wildlife and to prevent wind interference, the suitable sites should be more than 250m from woodland [35,36].  Resources Wales [42], SSSI in Scotland from Scottish Natural Heritage [43], airports from ShareGeo Open 220 [44], and existing wind turbines from the Virtual Wind Farm Model [32].
The total available areas were determined by taking the union of the suitable areas defined by constraints 3 and 4 and then intersecting with those of the other 8 constraints. Figure 2 shows the available area for wind turbines in each zone, which defines the land footprint constraints in the model. For comparison, Table 1 gives the available areas before and after application of the 10 constraints.

Temporal discretisation
The model needs to take into account simultaneously the long-term strategic decisions as well as short-term operational issues. The challenge in modelling is the need to represent energy storage using short time scales (i.e. not coarser than hourly intervals) to capture its dynamics. Considering a whole year planning horizon and formulating the model using contiguous hourly intervals result in a computationally intractable 230 model. Several methods to overcome this problem were discussed in [1]. The non-uniform hierarchical time discretisation and a decomposition method were applied in this paper to solve the model within an acceptable time.
Instead of representing time as contiguous hourly intervals, different time layers can be used: yearly intervals to model investment decisions, seasonal intervals to model seasonal variations (e.g. in demand and wind 235 availability), daily intervals to capture the difference between weekdays and weekends and hourly intervals for system balancing. This approach allows for a more efficient representation of time by exploiting the periodicity inherent in some of the system's properties. For example, instead of using 7 daily intervals per week, the periodicity in the demand data can be utilised, e.g. weekdays are likely to have similar demand profiles, which are likely to be different from those for weekends. Therefore, a particular day type can be 240 repeated a certain number of times and the demands can be represented as a sequence of repeated profiles, e.g. 5 repetitions of a weekday profile followed by 2 repetitions of a weekend profile can represent a week's worth of data. Since energy storage is one of the technologies being modelled, the storage inventories within and between time levels have to be linked, thus requiring additional variables and constraints. In this work, a single year was considered and a time discretisation of 4 different seasons in a year (spring, summer, autumn, 245 winter), 2 different day types in a week (weekday and weekend) and 24 hours in a day was used.

Spatio-temporal domestic road transport demand
The energy demand for transport in the domestic sector on a 1km square grid was estimated by Wang et al. [45], the data for which can be downloaded from [46]. They disaggregated the data for GB using the statistical data from the Living Costs and Food Survey [47] and census data. The data are "home-based" 250 rather than "road-based", i.e. the data were assigned to the homes of the drivers of the cars and vans rather than the stretch of road where the emissions were produced. This assumption does not affect the suitability of the data for this study because the 1km data were aggregated to the 16 transmission zones considered in this study. The 1km domestic transport data provided by Wang et al. were assumed to be demands for petrol. These were converted to demands for hydrogen by using the average fuel economies of petroleum 255 cars and fuel cell cars, the former being 49 mpgge and the latter being 79 mpgge [48]. Figure 3   The temporal distribution of demands was estimated from the statistical data set on traffic flows, which can be downloaded by month (TRA0305), by day of the week (TRA0306) and by time of day (TRA0307) 260 from the Department for Transport website [49]. These data were used to disaggregate the average hydrogen demand in each zone, shown in Figure 3(b), into 4 different seasons, 2 different day types (weekday and weekend) and 24 hours in a day. An example result of the disaggregation is shown Figure 4, which presents profiles show a distinct bimodal shape, corresponding to the morning and evening rush-hour periods, and the demands in the weekends peak around noon.

Spatio-temporal wind availability
The Department of Trade and Industry wind speed database [39] provides estimates of the annual wind speed throughout GB. Figure 5 shows the locations with annual average wind speed of at least 5m/s at 45m 270 above ground level.  In this study, the spatio-temporal wind speed data were obtained from the Virtual Wind Farm Model [32].
Historic weather data from 2014 were used to produce wind speeds for 8,760 hours for 10 different locations in each zone, which were then aggregated to determine the representative wind speed in each zone. In order to match the temporal discretisation used in the case studies, a representative profile for each day type and 275 each season is required. Averaging the profiles over all days of the same day type over all weeks in the same season resulted in the lose of variability in the data which may lead to an under-designed network. Therefore, for a conservative design, the daily profiles were chosen to be the most variable (defined as the one with the largest difference between the maximum and minimum wind speed) among all of the different days of the same day type over all weeks in the same season and, where possible, profiles for different day types were 280 chosen to be different to each other. Figure 6 shows the resulting wind speed profiles for zone 13; the profiles for other zones are not shown in the interest of space. Given the large uncertainty in the behaviour of the wind, the actual operation of the network will change to accommodate different scenarios as they occur but the network design must of course remain fixed. Using different wind profiles in the optimisation may result in different network designs, e.g.: flatter wind profiles, or ones that match more closely the demands, may 285 result in less installed storage capacity; the more the wind profiles are variable or incompatible with the demands, the higher the installed storage capacity is expected to be.

Network structure
As already mentioned, in this study GB is divided into 16 transmission zones. To illustrate the structure of the network, Figure 7 shows two example zones. In each zone, a number of wind turbines may be installed 290 in order to generate electricity. This can be converted to hydrogen, using electrolysers, which is used to fulfil transport demand. The hydrogen produced by the electrolysers is assumed to be at 20MPa, which is also the pressure at which it can be stored in underground caverns and pressurised vessels. Therefore, the hydrogen produced by electrolysers is labelled "Hydrogen at 20MPa" in Figure 7 and only this state can be stored in "Underground storage" and "Pressurised vessels". The hydrogen produced by the electrolysers or withdrawn 295 from storage must ultimately be delivered to customers via underground distribution pipelines to fuelling stations. The hydrogen supplied to customers may have been produced locally or imported into the zone from another zone within GB. Transport of hydrogen from one zone to another takes place in transmission pipelines, which based on the data provided by Yang and Ogden [50], have a maximum inlet pressure of 7MPa and it is also assumed that the distribution network requires hydrogen at this pressure. Therefore 300 a second hydrogen state, "Hydrogen at 7MPa", is produced when the "Hydrogen at 20MPa" is expanded using the "Expanders" technology, which also produces some electricity. This lower-pressure hydrogen can then be distributed to customers within the zone (it is assumed that the fuelling stations are equipped with compressors to dispense hydrogen at a pressure required by fuel cell cars) or transmitted to another zone (where the demand may be higher but with fewer wind turbines available). Any hydrogen received 305 from other zones may also be stored but this requires the hydrogen to be compressed from 7MPa back to 20MPa by using compressors, which also consume electricity. It was assumed that the pressure drop in the pipeline is negligible to restrict the number of hydrogen pressure levels to two, which reduces the size of the model; the validity of this assumption was confirmed when the results of the case studies were obtained (see Section 7.1, para. 3). Finally, the hydrogen may also be converted to electricity in fuel cells, where a 310 maximum inlet pressure of 7MPa is assumed. Energy may also be transmitted between zones in the form  storage or convert it to hydrogen and transmit the hydrogen to another zone; alternatively, the model may choose to reduce the excess electricity production by disengaging some or all of the wind turbines. Figure 7 is the Resource-Technology Network (RTN) diagram for the problem being modelled, which represents all of the possible energy pathways in the system. An RTN comprises 2 nodes: resources (usually drawn as circles) to represent any distinct material state, e.g. having a particular composition, temperature 320 and pressure, and technologies (usually drawn as rectangles) to represent processes that convert a set of input states to a different set of output states. More discussions about RTN can be found in [1,51]

Electrolysers
The electrolyser data were based on the report by the National Renewable Energy Laboratory (NREL) [55].
The maximum production rate is 50 tonnes H 2 per day (69.38MW hy ) and the efficiency is 50kWh electricity per kg H 2 (i.e. 67%). The unit capital cost is £31.56M and the annual O&M costs are 5% of the capital 355 cost. The operating pressure of the electrolyser was assumed to be 20MPa, therefore the generated hydrogen can be directly put into storage but needs to be expanded to 7MPa for transmission or distribution.

Fuel cells
This work considers a 41.63MW el solid oxide fuel cell (SOFC) with an efficiency of 60%. The cost of which was estimated from a 200kW el SOFC unit manufactured by Bloom Energy using a sizing exponent of 0.85.

360
The unit capital cost is £44.86M and the annual O&M costs as a fraction of capital cost per year is 6%.

Compressors
The electricity required to compress hydrogen from 7MPa to 20MPa was calculated to be 0.56kWh/kg H 2 by simulating a compressor train with interstage cooling in gPROMS ProcessBuilder [30]. Compressors of 7 different sizes were considered -the size of each one was determined based on the maximum injection rate of annual O&M costs were assumed to be 5% of the capital cost, which include changing oil regularly, replacing valves when needed among others. Table 2 summarises the characteristics of the compressors considered in the case studies; the storage IDs given in the third column are defined in Tables 6 and 7.

Expanders
Similar to compressors, expanders of seven different sizes were considered; the size of each one was based on the maximum withdrawal rate of each storage technology (discussed in Section 4.3). The electricity that can be recovered from the expansion of hydrogen from 20MPa to 7MPa was calculated to be 0.29kWh/kg H 2 by simulating a train of expanders with interstage heating in gPROMS ProcessBuilder [30]. The capital cost 375 was estimated from an expander with a rated output of 1MW el described in [56] using a sizing exponent of 0.80. The annual O&M costs were assumed to be 5% of the capital cost. The properties of the expanders used in the case studies are given in Table 3. study, since the demand density at the 1km level is available, it was assumed that the zones are connected by their centres of demand, obtained from the spatially-distributed demand at the 1km level and equations 1 and 2: where x z is the x-coordinate of the demand centre of zone z, x i is the x-coordinate of the centroid of 1km cell i, which is in zone z if i is in the set I z , andD i is the average demand in 1km cell i; the y-coordinate is 390 calculated similarly. The centres of demand are indicated by the cross-marks in Figure 3(a).
Since the underground storage facilities are not located at the centre of demand of each zone, the additional cost of transporting between the centre of demand and the underground storage was included in the cost of underground storage. The number of fuelling stations depend on the total demand in each zone whereas the length of the distribution network in each zone is estimated from the centres of demand and the distribution 395 of demand at the 1km level (discussed in Section 6).

Hydrogen pipeline
Yang and Ogden [50] provided the following properties of hydrogen transmission pipelines: maximum inlet pressure of 7MPa, outlet pressure of 3.55MPa and diameter of 100cm. Given these conditions, and the length and angle of inclination of the pipeline, it is possible to calculate the maximum flow of hydrogen 400 through the pipeline. The maximum steady-state flowrate occurs when the pressure difference across the pipeline is balanced by the frictional forces at the wall and the gravitational forces (for inclined pipelines).
The frictional forces depend on the length of the pipeline, the roughness of the wall and the square of the velocity of the fluid. Thus the force balance gives the velocity of the fluid in the pipeline and the maximum flowrate is obtained from F H2 = ρ H2 u H2 A pipe , where u H2 is the velocity of the fluid, ρ H2 is the density of 405 the fluid and A pipe is the cross-sectional area of the pipe.
There is an additional restriction on the throughput of the pipeline due to the velocity of the fluid, which should be lower than the erosional velocity. Typically, the pipeline is operated to ensure that the velocity is always lower than the erosional velocity by a specified margin (the erosional velocity margin). In this study, the velocity was restricted to be no more than 70% of the erosional velocity, calculated as: u er = 410 121.99/ ρ H2 .
Since the density of the hydrogen in the pipeline depends on the pressure, which drops from 7MPa at the inlet to 3.55MPa at the outlet, a more accurate approach to calculating the flowrate is to derive and solve a coupled set of partial differential and algebraic equations from differential mass, momentum and energy  Table 4. The shorter pipeline reaches its maximum flowrate (98kg/s) with an inlet pressure of just 4.73MPa, which indicates that the velocity in the pipe has reached 70% of the erosional velocity. Similar results can be seen for the longer, 230km, pipeline, where the simulation ended with the inlet pressure at 6.77MPa and a maximum flowrate of 82kg/s. Taking the lower flowrate of 82kg/s and multiplying by the lower heating value of hydrogen gives, as a conservative estimate, a maximum energy throughput of the transmission network 430 of 9.81GW hy . It may be noted that the erosional velocity is lower in the 230 km pipe than it is in the 48km pipe. This is due to the higher average pressure (since a much higher pressure drop is required to achieve the same flowrate as in the 48 km pipe), which results in a higher density and thus a lower erosional velocity. The capital cost of the pipeline comprises cost of materials, installation, right-of-way and miscellaneous.
Similar to the approach used in [50,61], the cost of materials was calculated as a function of the pipeline diameter, which is approximately £1.74M/km for pipelines with a diameter of 100cm. The pipelines were assumed to be installed next to the existing pipelines, hence the rights of way cost can be neglected. The installation cost is higher in urban locations than in rural places but in this study an average installation cost of £270k/km was used. Since the operating pressures of the electrolyser and hydrogen storage are much higher than the maximum pressure of the pipeline, it was assumed that no compression would be required 440 for transmission of hydrogen. A turbine at the pipeline inlet may be used to reduce the pressure to 7MPa and generate electricity. The annual O&M costs were assumed to be 2% of the capital cost.

Electricity transmission lines
Different high-voltage electricity transmission options were considered in the case study: alternating current overhead lines (HVAC OHL) and underground cables (HVAC UC) and direct current overhead lines (HVDC 445 OHL) and underground cables (HVDC UC). The properties of these technologies, which were estimated from [62,63], are summarised in  2. Injectability, which is the maximum rate at which gas can be injected into storage, and 3. Deliverability, which is the maximum rate at which gas can be withdrawn from storage.
In the case studies, overground and underground hydrogen storage facilities of different sizes were considered and their properties are described in the following sections.

Pressurised storage vessels
One of the advantages of this type of storage is its simplicity: the only requirement is a compressor and a pressure vessel. One of the disadvantages, however, is the low storage density, which depends on the storage 470 pressure. In general, higher pressure means higher capital and operating costs [64]. In the case studies, storage tanks of 3 different sizes, operating up to 20MPa, were considered, the characteristics of which are summarised in Table 6. Since compressors are considered as conversion technologies rather than components of storage technologies, the capital cost only includes the vessel and the cushion gas requirement. Cushion gas is the volume of gas required to be kept in the facility in order to maintain the operating pressure and cannot 475 be recovered until the facility stops operation. For compressed gas storage, the cushion gas requirement or "heel" was assumed to be 7.5% of the total capacity [65]. The annual O&M costs, which cover personnel and maintenance costs, were assumed to be 2% of the capital cost, which is the average of the values reported in [65]. Potential candidates for underground storage include salt caverns, depleted oil/gas fields and aquifers. There are a significant number of underground storage facilities, both at the operational and planning stage, in GB (cf. Figure 1 in the fact sheet provided by British Geological Survey [66]). Only 4 underground storage facilities, for which data are available in the literature, were considered in the case studies. Table 7 gives a summary of the characteristics of the 2 salt caverns (Aldborough and Warmingham) and 2 depleted oil/gas 485 fields (Humbly Grove and Rough) considered in the case studies. The data for the maximum available capacity, injectability and deliverability were obtained from the National Grid [67]. The capital costs, which were estimated from the report by Thistlethwaite et al. [68], include land/depleted field acquisition, cavern construction (leaching and brine disposal) for salt caverns, wells and above ground treatment facilities, connecting pipelines, cushion gas and so on. The cushion gas requirement was assumed to be 20% for 490 salt caverns and 45% for depleted oil/gas fields [65]. The annual O&M costs were assumed to be 2% of the capital cost [65]. Since the location of the underground storage does not coincide with the centre of demand of the zones where they are located, the additional cost for transporting the hydrogen between the underground storage and the centre of demand was also included in the costs. Other components that can further contribute to the costs but were not considered in the estimates include sub-surface analysis (e.g.

495
seismic data), control systems and planning and environmental approvals.

Model formulation
The model presented in this section is based on the general spatio-temporal energy systems model, STeMES, developed by Samsatli and Samsatli [1]. This section briefly describes the model and the adaptation of the general framework to include integrated wind-hydrogen-electricity networks. 500

Objective Function
The objective is to minimise the total annual costs comprising operating and capital costs of the network. . The annual costs of importing and exporting resources (e.g. buying and selling of electricity from the conventional grid) are given by equations 8 and 9, respectively, where V M r and V X r are the unit import and export costs of resource r, respectively (for exports, the cost is negative if the resource has value and can be sold); M rzhdt and X rzhdt are positive variables representing the rates of import and export, respectively, of resource r in zone z at hour h, day d and season t; n hd h is the duration of hourly interval h, n dw d is the number of times day d occurs in a week and n wt t is the number of repeated weeks in season s.
The total capital costs of the technologies can be determined from the product of the number of technologies and their respective unit capital cost. These are given by equations 10 to 13 for wind turbines, production 550 technologies, storage facilities and transport infrastructures, respectively, where C WT , C P p , C S s and C B b are the capital costs of a single wind turbine, production technology p, storage technology s and per unit distance of transmission infrastructure b, respectively (all £/technology installed except for C B b which is £/connection/km).

Resource balance
The operation of the network is governed by a resource balance given by constraint 14: where U rzhdt is the rate of utilisation of naturally-occurring energy resource r, M rzhdt is the rate of import of resource r from outside of the network into zone z (e.g. purchase of electricity from the conventional grid), P rzhdt is the net rate of production of resource r due to the operation of conversion technologies, S rzhdt 560 is the net rate of utilisation of resource r from storage in zone z, Q rzhdt is the net rate of transmission of resource r into zone z from other zones, D rzhdt is the demand for resource r (which, in the case studies, is zero for all resources except for hydrogen at 7MPa) and X rzhdt is the rate of export of resource r from zone z to outside of the network (e.g. selling of electricity to the conventional grid).

Resource availability 565
The electrical power that the turbine generators can extract from the wind depends on the efficiency of the wind turbines, which cannot be greater than the Betz Limit of 59.3% [69], the rotor sweep area and the wind speed. Equation 15 gives the wind power potential, in MW el , in zone z at hour h, day d and season t, u max Elec,zhdt = 0.5 × 10 −6 ηρ air π N EWT z R EWT,ave 2 + N NWT where η is the efficiency of the wind turbine (power coefficient), ρ air is the air density (taken to be 1.23kg/m 3 ), R EWT,ave is the average radius of the existing wind turbines, which was taken to be 35m from the data given 570 in [32], R NWT is the radius of the new wind turbines (50m) and v zhdt is the wind speed, which varies with both location and time (see Section 3.5). In the case studies, u max rzhdt is zero for all other resources.
The wind power potential, u max Elec,zhdt , is the upper bound on the rate of utilisation of electricity, U Elec,zhdt . That is, one cannot utilise more electricity than can be generated from the wind with the number of wind turbines installed in the cell. This is expressed more generally below.
Constraint 17 ensures that the number of existing turbines that are used as part of the network does not exceed the total number of existing turbines: where the parameter N EWT,tot z represents the total number of existing wind turbines in zone z.
The land footprint constraint guarantees that there is sufficient land area suitable for the new wind turbines assuming a spacing between turbines of 5 rotor diameters: where A max z is the total available area in zone z obtained by applying the 10 constraints discussed in Section 3.2.

Conversion technologies
The net rate of production of resource r, P rzhdt , is defined as: where P pzhdt is the rate of operation of technology p and α rp is the conversion factor, which is the net 585 production of resource r per unit operation of technology p (it is positive if r is produced and negative if r is consumed).
The production rate of conversion technology p in zone i is limited by the number of technologies of type p that are installed in zone i and the minimum and maximum capacities of a single technology: The number of technologies that can be established every year in the whole study region (i.e. build rate) can 590 also be constrained as follows: where BR p is the maximum build rate for technology p. Figure 9 shows the three stages involved in storing hydrogen: charging, maintaining and discharging, which are considered in the model as three storage tasks: "put", "hold" and "get". Similar to a conversion technology, 595 the efficiency, resource requirements and losses of each storage task can be specified through its conversion factor. It is also necessary to define the direction of the flow of resources, i.e. the source and the destination (abbreviated in the formulation as "src" and "dst", respectively). The "put" task transfers the hydrogen from zone z (source) to the store (destination). The "hold" task maintains the hydrogen in the store; if there are any losses while it is being maintained (e.g. hydrogen gas escaping to the atmosphere or in the case of 600 electricity storage, the losses are in the form of heat) then the source is the store and the destination is the zone. The "get" task retrieves the hydrogen from storage (source) and delivers it to zone z (destination).

Storage technologies
In this work, compressors and expanders are explicitly modelled as conversion technologies. Therefore it is not necessary to define the electricity requirement (or electricity generation in the case of expanders) of the "put" and "get" tasks. ToTtransmission/ distribution Figure 9: Three stages ("tasks") for storing hydrogen: "put", "hold" and "get".
The net rate of utilisation of resource r from storage is given by the rates of operation of all of the "put", "hold" and "get" tasks multiplied by their conversion factors for resource r: S rzhdt = s S put szhdt σ put sr,src + S hold szhdt σ hold sr,dst + S get szhdt σ get sr,dst The rates of charging and discharging the store cannot exceed the injectability, s put,max s , and deliverability, s get,max s , respectively. The inventory balance for the store also depends on the rates of the three tasks, but this time multiplied by 615 the conversion factor for the opposite flow direction: The rate of operation of the "hold" task is the inventory level from the previous time interval, divided by the length of the time interval: S hold sz,1,dt = I 0,sim szdt /n hd 1 ∀s ∈ S, z ∈ Z, d ∈ D, t ∈ T (26) Here, I 0,sim szdt is the initial inventory level for the start of the "simulated cycle" for day d in season t. It is calculated (equation 29) so that the inventory levels in the simulated cycle will correspond to the average 620 inventory levels over all occurrences of day d and all weeks in season t, so that costs and resource requirements, which depend on the inventory levels, are calculated correctly. See the appendix in [1] for a detailed derivation.
Equations 25 and 27 can be rearranged to give the inventory balance in a more familiar form: The left hand side is the rate of change of the inventory level, the first and third terms on the right hand side are rates of addition and withdrawal from storage and the second term on the right hand side is the rate of loss of resource from storage.
The variable I 0,sim szdt in equation 26 represents the initial inventory at the start of the "simulated cycle" for day d in season t. In order to make the model more efficient, only one of the identical days of each day type is considered and only one of the identical weeks of each season is considered. As the costs and resource 630 requirements for storage depend on the inventory levels, the day and week that is included in the optimisation, the "simulated cycle", should correspond to a day with inventory levels equal to the average inventory levels over all of the day types and all of the weeks in the particular day type, d, and season, t. If I 0,act szdt is the initial inventory for the first occurrence of day type d and the first week of season t, then the initial inventory for the simulated cycle is given by: δ d szdt is the change in inventory over one day, for day type d in season t and n dw d is the number of identical days in day type d. Therefore adding integer multiples of δ d szdt to I 0,act szdt gives the initial inventory on subsequent days (of the same day type) and adding (n dw d − 1)/2 multiples of δ d szdt shifts the initial inventory to the average level over all occurrences of that day type. δ t szt and n wt t are defined as the change in inventory over a week and the number of weeks in season t; they are combined similarly to shift the inventory forward to 640 the central week in season t. These relationships are derived in more detail in the appendix of [1]. δ d szdt and δ t szt are defined below.
The change in inventory over the whole year can also be calculated and used to ensure that there is no accumulation of inventory over a year: The initial inventories for each day type are also related: adding n dw d multiples of δ d szdt to the initial inventory 645 for day type d gives the initial inventory of day type d + 1. The equations below link the initial inventories for subsequent day types and seasons.
The following set of constraints ensures that the inventory does not exceed the maximum capacity of storage, s hold,max s . For computational efficiency, the periodicity in system properties can be exploited, i.e. the inventory increases or decreases by the same amount, δ d szdt , each repeated day in day type d in season t 650 and by the same amount, δ t szt , each repeated week in season t. Therefore, instead of writing the constraints for every contiguous hour, it is sufficient to write the constraints only for the first and last instance of each repeated day type and the first and last week of each season, as given by constraint 36.
Note that the above is a short hand for the four constraints formed by using either a positive or negative sign for each of the ± symbols.

655
In the case studies, it was assumed that each storage should be matched with a dedicated compressor and expander, each of which is specifically sized for each storage (see Tables 2 and 3). Therefore, two constraints were added for each row in Tables 2 and 3. Those for the first row (i.e. for the large H 2 storage tank) are shown below.

Transport technologies 660
The net rate of transport of resource r into zone z from other zones, Q rzhdt , is the difference between the incoming and outgoing flow rates given by the first and second terms on the right hand side of equation 39, respectively: Q lzz hdt is the rate of operation of transport mode l from zone z to zone z during hour h, of day type d in season t. The conversion factors,τ lrf , are the net flow of resource r into the source (f = src) and destination (f = dst) zones per unit rate of operation of the transport mode: they are negative if the flow is out of the zone and positive if the flow is into the zone.τ lrf are defined similarly but are also per unit distance between the zone, hence they are multiplied by the distance between the zones, where x z and y z are the x-and y-coordinates of the centre of demand of zone z, and the rate of operation of the transport mode; they are mainly used to represent distance-dependent losses. By convention,τ lr, dst 670 is set to 1 for the resource being transported so that Q lzz hdt is the rate of resource r entering zone z from zone z (and q max l is the maximum rate of transport of the resource, see equation 40 below). For a lossless transport mode, with no resource requirements,τ lr, src is set to −1 for the transported resource and all other elements are set to zero.
The rate of operation of transport mode l, Q lzz hdt , is limited by the maximum rate, q max l , and the number 675 of infrastructures of type b, N B bzz , established between zones z and z : Here, LB lb is a binary parameter with a value of 1 if transport mode l is supported by infrastructure b, 0 otherwise. For computational efficiency, the infrastructure links can be only built between adjacent (or neighbour) zones. The transport of resource r over long distances can be achieved by making several neighbour-to-neighbour transfers along the route between the source and destination zones. The binary 680 parameter, ν zz , is used to indicate neighbour zones: ν zz = 1 if zone z is adjacent to zone z , 0 otherwise.
The total flow rate of all resources being transported along an infrastructure is also constrained by its Finally, the following constraint applies for bidirectional transmission lines, such as the pipelines and cables considered in this study:

Import and export
The rates of import and export can be constrained by specifying the maximum rates of import and export, m max rz and χ max rz , respectively.

Distribution network
It is assumed that customers will purchase hydrogen from a number of fuelling stations distributed throughout 690 each zone, much as they currently do for petrol and diesel. Indeed, it would make sense for some existing petrol stations to convert to hydrogen, in which case the distribution network can be designed in detail.
In this study, the exact locations of the fuelling stations have not been considered so the properties of the distribution network must be approximated from the spatial distribution of demands. In the earlier study of Almansoori and Shah [57,58], the fuelling stations were assumed to be uniformly distributed within each which they assumed to be half the length of the square [57]. These assumptions were not appropriate for the following reasons: first, Figure 3(a) shows that the demand for hydrogen is far from uniform; second, it can be shown that the average distance from the centre of a square to all points in the square is in fact where L sq is the length of the square. Figure 10(a) compares the effect that the assumption of uniform demand has on the estimated length of the distribution network: generally, 705 the length is overestimated and in one zone the error is roughly 90%. In this study, it was assumed that the hydrogen is distributed by pipeline, so the total distribution network distance in each zone can therefore be estimated from the 1km demand data using equation 45: where L network z is the length of the distribution network in km,D i is the demand density at the 1km level and F S is the capacity of a single fuelling station, which was assumed to be 1,500 kg/day (2.08MW hy ) and 710 is the size of a facility described in the report by NREL [70].
The total network distance is required to cost the distribution network and the demand-weighted average distance from the centre of demand to all of the points in the cell (L network z /F S) can be used to estimate distribution losses in cases where the losses are proportional to distance from the supply point to the customer (one can easily perform more general calculations to arrive at a representative number to use in the loss 715 constraints).
The number of fuelling stations in each zone, N FS z , can be determined by dividing the total demand in each zone by the capacity of a single fuelling station and then rounding up to the nearest integer (equation 46), the results of which are shown in Figure 10

Sizing
The pipeline diameter is initially assumed to be 20cm but it is important to ensure that the capacity of the distribution network exceeds the peak demand in each zone. A detailed design of the distribution networks is beyond the scope of this study but it is possible to perform a similar analysis as was given in section 4. 10, which are at 69%, 144% and 96% of the erosional velocity and therefore require larger-diameter pipes to connect the demand clusters to the centre of demand. Increasing the pipe diameter to 30 cm brings these numbers down to 26%, 44% and 34% respectively. For additional robustness, the pipes in zone 8 were sized at 35 cm in order to bring the velocity down below half of the 70% limit, in line with all of the other zones, which can support a doubling of the hydrogen demand.

Costing
The total annualised cost of the distribution network is estimated from the network length and the number of fuelling stations determined by equations 45 and 46, as well as the diameter of the pipeline as determined above and a standard capital charge factor of 3. The unit cost of the pipe was estimated to be £348k/km, £437k/km and £498k/km for 20, 30 and 35cm diameter pipes using a similar approach to that of Yang and 755 Ogden [50] and Parker [61] and the unit cost of the fuelling station was approximated to be £3.03M based on the values reported in an NREL report [70].This yields a total annualised cost for the distribution network of £17.1bn/yr. Although the cost may be reduced by optimising the distribution network, due to the very large size of the model that will result if the distribution is to be included in the optimisation problem, which is unlikely to be solvable using available computing resources, only networks comprising generation/conversion, 760 storage and transmission technologies are optimised and the estimated cost for the distribution is simply added to the total cost.

Results and discussion
A number of different scenarios were considered, using the data described above, in order to investigate the different optimal configurations of the network. The base case includes all of the technologies but does not 765 use existing wind turbines. In subsequent runs, the value of various network components was determined by excluding them from the optimisations and comparing the resulting cost to the base case. Another scenario also allows the existing wind turbines to be used. Although it is straightforward to consider imports and exports, as discussed in Section 5.2.6, none of the scenarios allowed them.
The overall conclusion is that both storage and transmission are required to obtain a feasible solution, i.e.

770
without either of these, the energy demands cannot be met at all times. Transmission is needed because some zones in the south of England do not contain sufficient suitable land area for wind turbines in order to meet the demands independently, even with storage, i.e. energy must be supplied by other zones to meet the demands. Storage is required even if all zones are connected by unlimited transmission capacity because the intermittency of the wind availability is too high to meet the total demand at all times.

775
Of course, these results are to be expected but the interesting questions that the model is designed to answer are: how much storage is required and where should it be located; what type of transmission is required and which zones should be connected? Naturally, in addition to these, the model determines the number and location of all of the other technologies, such as electrolysers, fuel cells etc.
The different cases examined are described in the subsequent sections.

Base case
For the base case, all of the technologies in the database were considered but without the availability of the existing wind turbines, assuming that their generation capacity is already allocated to satisfying other demands. The optimal structure of the resulting network is shown in Figure 11(a). The symbols, with numbers, in each zone represent the type and number of technologies in the zone and the lines represent The cost breakdown is shown in Figure 11 To these costs should be added the distribution network costs of £17.1bn/yr, as estimated in Section 6.2.
However, since the distribution costs are the same in all cases, they are not included in the comparisons made in the following subsections.

800
The model also determines the hourly operation of each technology. For example, Figure 12 shows the operation of the pipeline transmission network at different times during weekdays in summer. At times of high demand, zone 15, where the Humbly Grove underground storage is located, supplies a large amount of hydrogen to other zones that are connected to the transmission pipeline. Zones 7 and 16, which have significant generating capacity, supply hydrogen to other zones through the transmission pipeline at some 805 times and also receive hydrogen from other zones through the pipeline at other times. Zones 9 to 13 satisfy all of their domestic transport demand through the pipeline at all times. The maximum flow of hydrogen through the pipeline at peak time is 2.85 GW hy and the pressure drop was calculated using gCCS to be 0.06 MPa. Therefore, the decision not to model the booster compressors in the pipeline was justified.  During times of low demand, the excess production of hydrogen is being stored in the underground storage;

810
an example of such instance is presented in Figure 13, which shows the operation of the transmission pipeline at 01:00h during weekdays in spring. The hydrogen being generated in zones 7 and 16 is being transmitted through the pipeline to zone 15. Because the transmitted hydrogen is at a lower pressure of 7MPa, a compressor is required to raise it up to the storage pressure of 20MPa. As can be seen in Figure 13, the Humbly Grove underground storage in zone 15 is equipped with both a compressor and an expander, 815 the former is needed when the hydrogen from transmission is stored while the latter is required when the hydrogen from storage is transmitted to other zones or distributed within the zone to meet demands.

Computational statistics
The model was implemented in AIMMS 3.12 and the decomposition method described in [1] was employed  and designing the transport infrastructure with a fixed set of conversion and storage technologies. The decomposition algorithm terminates when either stage 2 or 3 no longer improves the objective function. Table   8 shows the computational statistics for each iteration of the decomposition procedure. As the objective function does not improve in stage 3, the decomposition algorithm terminates at the third iteration and the 850 solution for stage 2 is taken to be the optimal solution.
As can be seen from Table 8 In case 2, the existing wind turbines can be used as part of the network. As can be seen in Figure 15   cost of the existing wind turbines does not contribute towards the total cost; only their O&M costs do), the conversion (electrolysers, expanders and compressors) and transmission costs are similar (i.e. within the 3% relative tolerance used as the terminating criterion when solving each stage of the model) and the cost of storage is 20% higher. The cost of this network is 7% lower than that of the base case.

The value of underground storage
This case is the same as the base case except that the use of underground storage is not allowed. Figure   885 15(c) shows the resulting optimal network structure. completely relies on other zones for the satisfaction of its demand. In contrast, zone 11, which does not have any technology in the base case, now has a significant capacity for generation and storage and also supplies hydrogen to its neighbouring zones through the transmission pipeline. As can be seen in Figure 15(d), 895 without underground storage the cost of the network is 25% higher than that of the base case: the wind turbines, conversion, storage and transmission are more expensive by 29%, 23%, 32% and 20%, respectively.

The value of pipeline transmission
In this case, the effect of not being able to transmit hydrogen through pipelines was examined. Figure 15( to other zones. In the base case, zones 9 to 13 do not require any conversion or storage capacity because they satisfy their demands using hydrogen from the transmission pipeline. Without the hydrogen pipeline, these zones now need electrolysers to generate hydrogen using either the electricity from the local wind turbines (e.g. zone 11) or the overhead electricity cables.
Figure 15(f) shows the cost breakdown for this network. Compared to the base case, the cost of transmission 915 is lower by 31% but the cost of wind turbines, conversion technologies and storage facilities are higher by 7%, 17% and 57%, respectively. Overall, this network is 11% more expensive than that of the base case.

The cost of underground electricity cables
Overhead electricity lines have a negative visual impact on the landscape, therefore in this case, the additional cost of using underground electricity cables as the the only option for transmission was determined.

920
The resulting optimal network structure is presented in Figure 15(g). All zones have electrolysers and storage facilities that are equipped with expanders. Except for zone 9, which obtains its electricity from the trans-

Conclusions
This paper presented an MILP model for the optimal design and operation of integrated wind-hydrogenelectricity networks. The general modelling framework, STeMES, was extended to include a footprint con-935 straint for wind turbines and a more detailed representation of the available wind energy. The area available for siting wind turbines was determined, using GIS, by applying 10 environmental and technical constraints, such as buffer distances from urban areas, rivers, roads, airports, woodland and so on. If the use of the existing wind turbines is permitted, the total cost of the network is 7% cheaper than the base case, despite the cost of wind turbines being 22% lower. This indicates that not all of the existing wind farms are in the ideal location and require more expensive storage in order to utilise them in the energy network. Without any underground storage, the optimal network is 25% more expensive than the base case. This indicates how important large, seasonal storage facilities are. Without hydrogen transmission pipelines, the optimal network, which uses HVAC overhead lines, is 11% more expensive than the base case. The transmission is less expensive but fuel cells and more storage are required, more than offsetting the savings in the transmission network cost. Lastly, if transmission is restricted to underground electricity cables, the cost of the network is 37% higher than the base case.
The results of the case studies are naturally dependent on the quality of the input data used, especially the 970 costs of the technologies, which were assumed to be independent of location but could vary widely across Great Britain due to a number of factors such as installation costs. It is straightforward to include these dependencies in the model but the availability of data is the main difficulty. While every effort was made to obtain reliable data, there is still significant uncertainty in their values. Therefore, a useful future step would be to perform sensitivity analyses. Moreover, the results are of course limited by the technologies that 975 are currently in the database and they may change with the inclusion of more technologies, such as batteries and electric vehicles -these are all planned future extensions to the model. In addition, the model will be extended to include other networks such as natural gas and heat along with the demands from other sectors.
Pipeline storage and hydrogen injection into the natural gas grid will also be investigated.